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1  ABSTRACT 


1  Abstract 

SERDP  project  MR- 1664  entitled  “Isolating  and  Discriminating  Overlapping  Signatures  in  Clut¬ 
tered  Environments”  is  approximately  halfway  complete.  Significant  progress  has  been  made  in 
working  toward  the  original  objectives  of  the  project.  Three  new  methods  for  localizing  multiple 
sources  in  close  proximity  using  EMI  data  have  been  developed  and  tested.  Specifically,  these 
methods  are: 

•  A  multiple  dipole  search  method  based  on  a  gradient  search  algorithm  utilizing  an  analytical 
Jacobian  (see  Sec.  4.2) 

•  A  combined  Joint  Diagonalization  (JD)  and  Orthonormalized  Volume  Magnetic  Source  (ON- 
VMS)  method  (see  Sec.  4.3) 

•  Source  localization  based  MUSIC  algorithm  applied  to  EMI  data  (see  [1]) 

Canonical  targets  of  various  shapes,  sizes  and  material  parameters  have  been  fabricated  (see 
Sec.  4.1.  Data  acquired  from  these  targets  as  well  as  standard  UXO  targets  has  been  acquired 
by  the  TEMTADS  and  MPV2  instruments  in  many  multitarget  configurations  (see  Sec.  4.1).  The 
methods  developed  to  date  under  this  MR- 1664  have  been  able  to  isolate  and  discriminate  up  to 
six  targets  simultaneously  in  the  case  of  lab  data  (see  Sec.  5).  The  JD  method  is  able  to  almost 
instantaneously  provide  a  good  estimate  for  the  number  of  distinct  targets  in  the  EMI  data.  After 
this  estimate  is  obtained,  the  first  or  second  methods  delineated  above  (and  described  below)  are 
used  to  invert  for  the  parameters  of  the  N  identified  targets. 
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2  OBJECTIVES 


2  Objectives 

The  objectives  of  the  proposed  work  are  to  mitigate  the  effects  of  discrete  metallic 
clutter  from  EMI  data,  and  the  simultaneous  discrimination  of  multiple  anoma¬ 
lies  in  close  proximity  with  overlapping  EMI  responses.  Specifically,  we  aim  to 
develop  a  robust  approach  for  processing  high  quality  data  from  next  generation 
EMI  instruments  to  mitigate  the  effects  of  clutter  by  isolating  their  EMI  signatures 
after  locating  them  individually,  and  then  discriminate  UXO  from  non-UXO  tar¬ 
gets  at  highly  contaminated  sites  that  include  challenging  terrain,  vegetation,  and 
geology  using  rigorous  models  that  may  include  interaction  effects. 

1.  Develop  an  TV-target  estimator,  able  to  provide  estimates  of  the  number  of 
targets  present  in  the  sensors  held  of  view  along  with  their  locations  and 
orientations,  without  resorting  to  computationally  expensive  optimizations. 
This  locator  will  utilize  high  quality  EMI  data  without  solving  the  traditional 
ill-posed  inverse  scattering  problem. 

Here,  two  approaches  for  this  TV-target  locator  will  be  reported  (Sec.  4): 

(a)  A  multiple  dipole  search  method  based  on  a  gradient  search  algorithm 
utilizing  an  analytical  Jacobian  (see  Sec.  4.2) 

(b)  A  combined  Joint  Diagonalization  (JD)  and  Orthonormalized  Volume 
Magnetic  Source  (ONVMS)  method  (see  Sec.  4.3) 

(c)  Source  localization  based  MUSIC  algorithm  applied  to  EMI  data  (not 
reported  here,  see  [1]) 

2.  Formulate  robust  classifiers  that  segregate  TV  targets  into  UXO  and  non-UXO. 
This  classifier  will  consider  model  based  parameters  such  as  volume  esti¬ 
mates  and  HAP  estimate  of  polarizability  and  may  be  based  on  rapid  statisti¬ 
cal  techniques  such  as  Support  Vector  Machines  (SVM). 

Significant  progress  on  this  has  been  made.  For  example,  our  group  recently 
achieved  a  very  encouraging  score  (no  false  negatives  and  very  few  false 
positives)  on  our  dig  list  for  the  Camp  Butner  demonstration  site  using  an 
Expectation  Maximization  (EM)  algorithm  (see  [2]). 

3.  Discriminate  UXO-like  targets  using  rigorous  (NSMS,  SEA)  models  which 
explicitly  include  coupling  between  targets  if  required  (see  Sec.  4.2-Sec.  4.3). 
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Unmatched  high  quality  data  from  advanced  instruments  in  the  time  domain 
will  be  used  (Sec.  4.1).  These  instruments,  specifically  in  this  report  the  TEM- 
TADS  and  MPV2  instruments,  provide  tremendous  data  diversity  including  vector 
data,  monostatic  and  multistatic  measurements,  allowing  us  to  assess  the  benefits 
of  data  diversity.  As  well,  rigorous,  physically  complete  models  (see  Sec.  4.2- 
Sec.  4.3)  are  being  developed  to  discriminate  targets  and  anomalies  based  on  their 
extracted  EMI  model  parameters.  These  high  fidelity  models  reduce  to  dipole 
models  but  can  also  provide  much  more  information  than  traditional  dipole  mod¬ 
els  alone.  The  JD  method,  for  example,  trades  knowledge  of  the  target’s  position 
and  what  the  amplitude  of  its  polarizabilities  are  for  a  rapid  determination  of  the 
shape  of  those  polarizabilities. 


3  Background 


EMI  data  acquired  over  subsurface  anomalies  must  be  processed  and  inverted  in 
order  to  discriminate  innocuous  scrap  metal  or  geology  from  munitions  and  ex¬ 
plosives  of  concern  (MEC).  At  highly  contaminated  sites,  processing  data  con¬ 
taining  responses  from  disperse  metallic  clutter  and  overlapping  anomalies  is  a 
difficult  task  both  theoretically  and  computationally.  Response  from  shallow  yet 
small  metallic  clutter  can  often  obscure  deeper  targets  of  interest  (for  example, 
see  Figs.  28-30),  and  the  electromagnetic  induction  (EMI)  response  from  targets 
in  close  proximity  may  affect  each  other  due  to  mutual  coupling  [3,  4]  and  con¬ 
fuse  or  confound  inversion  schemes.  These  realities  often  lead  to  misclassification 
at  best,  or  a  complete  failure  of  the  routine  at  worst.  Both  outcomes  raise  the 
probability  of  false  alarm  and  the  probability  of  a  missed  MEC. 

Prior  approaches  aimed  at  addressing  the  topic  of  discriminating  multiple  over¬ 
lapping  signals  have  serious  limitations.  Some  studies  cast  the  problem  as  a  blind 
source  separation  (BSS)  problem  [5-7].  These  studies  assumed  no  interaction  be¬ 
tween  anomalies  which  interaction  is  known  to  exist  [8],  were  constrained  to  a 
maximum  of  two  targets  [5],  were  simulations  that  had  prior  knowledge  of  the  tar¬ 
get  [6],  or  a  combination  of  similar  limitations.  Other  studies  used  data  from  mul¬ 
tiple  instruments  (such  as  simultaneous  EMI  and  GPR)  which  can  be  costly  and  is 
not  generally  available  [9-12].  ESTCP  project  UX-0415  attempted  to  use  a  pair  of 
dipoles  to  accommodate  2  targets  and  an  iterative  approach  for  isolating  multiple 
dipole  sources.  However  other  work  [13]  has  shown  that  even  single  multi-part 
UXO’s  require  multiple  dipoles  for  effective  representation,  and  an  optimization 
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may  not  be  able  to  distinguish  a  single  heterogeneous  UXO  from  multiple  UXOs. 

A  recent  SERDP  project,  MM- 1637,  Selecting  Optimal  Models  for  Inverting 
EMI  Data  concentrates  on  model  selection  and  data  quality  to  determine  whether  a 
given  data  set  should  be  modeled  as  one  or  two  targets.  They  then  use  a  probability 
based  inversion  scheme  to  invert  the  data  based  on  their  assumed  number  of  tar¬ 
gets.  Our  approach  differs  from  these  approaches  in  fundamental  ways.  First,  we 
are  using  analytical  techniques  to  mitigate  the  noise  due  to  disperse  clutter  which 
problem  MM- 1637  does  not  address.  Second,  we  are  developing  methods  to  locate 
N  targets  simultaneously  as  a  initial  step  either  prior  to  or  during  inversion  (Sec.  4) 
whereas  MM- 1637  only  addresses  the  cases  of  one  or  two  objects.  This  TV-target 
locator  (the  JD  method  (see  Sec.  4.3)  and  MUSIC  target  localizer  [1]  avoid  time 
consuming  and  computationally  complex  inversion  stages  that  MM- 1637  does  not 
avoid.  Finally,  we  have  access  to  the  latest  generation,  well  located  instruments 
which  provide  high  quality  data  from  realistic  held  conditions. 

In  these  studies,  the  objects  are  generally  modeled  using  simple  dipole  mod¬ 
els  to  accommodate  two  targets  and  the  targets  are  isolated  via  a  laborious  time 
consuming  inversion  approach.  Problems  during  the  inversion  phase  also  arise 
when  the  digital  geophysical  data  are  very  low  quality  and  simple  phenomenolog¬ 
ical  models  fail  to  accurately  represent  overlapping  EMI  signals  from  anomalies 
in  close  proximity.  In  addition,  many  single  heterogeneous  targets  may  be  well 
modeled  by  two  dipoles,  further  complicating  inversion  processes  based  solely  on 
dipoles.  In  these  cases,  the  parameters  of  the  phenomenological  model  (which 
are  used  as  discrimination  features  in  statistical  decisions  algorithms)  become  un- 
separable  and  as  a  result  the  classification  is  unreliable.  Our  approach  detailed  in 
Sec.  4.2  uses  multiple  dipoles  seeded  into  the  subsurface  space,  but  uses  the  ana¬ 
lytical  Jacobian  to  simultaneously  speed  up  the  gradient  search  and  provide  higher 
confidence  in  the  results. 

As  an  example  of  EMI  data  with  overlapping  responses,  Fig.  26  in  Appendix  A.2 
shows  typical  frequency  domain  data  from  a  GEM-3  acquired  over  multiple  buried 
targets  [14].  Figure  30  (top  left  subfigure)  also  show  GEM-3  data  from  a  105mm 
UXO  obfuscated  by  a  layer  of  discrete  clutter  affixed  to  plywood.  Even  though 
we  do  not  plan  to  use  the  GEM-3  for  this  research,  evident  in  this  data  is  the  need 
for  robust  processing  which  identifies  and  classifies  discrete  anomalies  and  higher 
quality  data  addressing  both  positional  and  instrument  noise. 

The  challenges  involved  with  the  discrimination  of  anomalies  in  close  proxim¬ 
ity  are  significant  and  persistent.  In  all  prior  cases,  data  quality  in  terms  of  the 
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signal  to  noise  ratio  (SNR),  positional  inaccuracies,  and  data  completeness  (e.g. 
survey  spacing),  has  been  a  significant  factor  which  limits  the  success  of  the  vari¬ 
ous  models. 

4  Technical  Approach 

4.1  Canonical  Targets 

This  project  shares  one  early  goal  with  two  other  projects  (MR- 1662  under  SAIC, 
and  MR- 1669  under  BAE)  which  started  in  the  same  FY  as  this  MR- 1664:  to 
fabricate  a  set  of  canonical  targets  as  a  basis  and  baseline  for  measurements  and 
testing.  These  targets  have  been  fabricated  and  consist  of  two  sets  of  spheroids 
and  spheres  combined  together  under  the  stewardship  of  Dan  Steinhurst  at  NOVA 
Research:  an  old  set  made  by  Dr.  Ben  Barrowes  while  at  MIT  in  2003  and  consists 
of  smaller  spheroidal  and  spherical  targets  under  25cm  in  the  longest  dimension 
(see  Table  1  and  Fig.  1),  and  a  new  set  fabricated  in  2009  and  2010  for  these 
projects  subcontracted  through  Delpsi,  FFC  (see  Table  2).  Figures  2-3  shows  the 
spheroids  from  the  second,  larger  set.  Details  on  these  targets  are  in  Sec.  4.1.  Data 
collected  from  the  TEMTADS  instrument  over  these  canonical  targets  is  summa¬ 
rized  in  Sec.  5.1  while  data  collected  over  multitarget  UXO  configurations  with 
the  MPV2  instrument  in  included  in  Sec.  5.3. 
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comp. 

type 

axis  (2a) 

axis  (2b) 

nominal 

e  =  b/a 

SI 

Iron 

sphere 

90.62mm 

90.62mm 

1 

S2 

Steel 

PS 

30.02mm 

182.19mm 

6 

S3 

Steel 

PS 

29.94mm 

90.28mm 

3 

S4 

Steel 

PS 

14.97mm 

90.76mm 

6 

S5 

Steel 

OS 

29.32mm 

4.56mm 

1/6 

S6 

Steel 

OS 

29.59mm 

9.65mm 

1/3 

S7 

Steel 

sphere 

29.87mm 

29.87nmi 

1 

S8 

Steel 

OS 

89.85mm 

28.39mm 

1/3 

S9 

Steel 

OS 

89.95mm 

15.32mm 

1/6 

A1 

A1 

PS 

30.17mm 

1 80.23mm 

6 

A2 

A1 

PS 

29.9mm 

91.29mm 

3 

A3 

A1 

PS 

15.04mm 

91.14mm 

6 

A4 

A1 

OS 

29.36mm 

4.34mm 

1/6 

A5 

A1 

OS 

29.36mm 

8.88mm 

1/3 

A6 

A1 

sphere 

29.91mm 

29.91mm 

1 

A7 

A1 

OS 

89.92mm 

30.38mm 

1/3 

A8 

A1 

OS 

89.98mm 

15.94mm 

1/6 

Table  1:  Physical  dimensions  of  steel  and  aluminum  spheroids.  PS  =>■  Prolate  Spheroid. 
OS  =>  Oblate  Spheroid. 


Figure  1:  Spheroid  collection  used  for  testing.  Specifications  are  listed  in  Table  1. 
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comp. 

type 

nominal 
axis  (2 a,  cm) 

nominal 
axis  (2b,  cm) 

nominal 

e  =  b/a 

SSI 

steel 

sphere 

10 

10 

1 

SS2 

steel 

sphere 

15 

15 

1 

SP1 

steel 

sphere 

10 

20 

2 

SP2 

steel 

prolate 

6.67 

20 

3 

SP3 

steel 

prolate 

5 

20 

4 

SP4 

steel 

prolate 

4 

20 

5 

SP5 

steel 

prolate 

6 

30 

5 

SP6 

steel 

prolate 

12.5 

50 

4 

SP7 

steel 

prolate 

10 

50 

5 

SP8 

steel 

prolate 

10 

30 

3 

SP9 

steel 

prolate 

15 

30 

2 

SP10 

steel 

prolate 

7.5 

30 

4 

SOI 

steel 

oblate 

20 

8 

2/5 

S02 

steel 

oblate 

20 

4 

1/5 

S03 

steel 

oblate 

20 

2.67 

1/7.5 

S04 

steel 

oblate 

20 

2 

1/10 

AS1 

aluminum 

sphere 

10 

10 

1 

AS2 

aluminum 

sphere 

15 

15 

1 

API 

aluminum 

sphere 

10 

20 

2 

AP2 

aluminum 

prolate 

6.67 

20 

3 

AP3 

aluminum 

prolate 

5 

20 

4 

AP4 

aluminum 

prolate 

4 

20 

5 

AP5 

aluminum 

prolate 

6 

30 

5 

AP6 

aluminum 

prolate 

12.5 

50 

4 

AP7 

aluminum 

prolate 

10 

50 

5 

AP8 

aluminum 

prolate 

10 

30 

3 

AP9 

aluminum 

prolate 

15 

30 

2 

AP10 

aluminum 

prolate 

7.5 

30 

4 

AOl 

aluminum 

oblate 

20 

8 

2/5 

A02 

aluminum 

oblate 

20 

4 

1/5 

A03 

aluminum 

oblate 

20 

2.67 

1/7.5 

A04 

aluminum 

oblate 

20 

2 

1/10 

Table  2:  Physical  dimensions  of  larger  set  of  steel  and  aluminum  spheroids  fabricated  in  2009- 
2010.  First  letter  of  designation  is  material  (A=Aluminum,  S=steel),  second  letter  is  shape 
(S=Sphere,  P=Prolate,  0=Oblate),  and  a  number  designation. 
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(a)  SSI  and  AS  1 


(c)  SP1  and  API 


Prolate  Spheroid 
20  cm  x  5  cm  x  5  cm 


(e)  SP3 


(g)  SP5  and  AP5 


•o 

(b)  SS2  and  AS 


(h)  SP6 


Figure  2:  Larger  steel  spheroids  SS1-SP6.  Specifications  are  listed  in  Table  2. 
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i 


(a)  SP7 


(c)  SP9 


(e)  SOI  and  AOl 


(g)  S03 


(b)  SP8 


(d)  SP10 


(f)  S02  and  A02 


(h)  S04 


Figure  3:  Larger  steel  spheroids  SP7-S04.  Specifications  are  listed  in  Table  2. 
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4.2  Multiple  Dipole  Inversion  Using  Jacobian  Gradient  Search 

The  simultaneous  detection  and  identification  of  multiple  targets  using  electro¬ 
magnetic  induction  time-domain  sensors  remains  a  challenge  due  to  the  fast  decay 
of  the  magnetic  field  with  sensor-target  distance.  For  example,  the  signal  from  a 
weak  yet  shallow  target  or  clutter  item  can  overshadow  that  from  a  much  larger  yet 
deeper  UXO,  potentially  resulting  in  erroneous  localization  and/or  identification. 
We  propose  in  this  paper  a  method  based  on  the  Gauss-Newton  algorithm  for  the 
inversion  of  multiple  targets  within  the  field  of  view  of  sensors  operating  at  electro¬ 
magnetic  induction  -  EMI  -  frequencies  (tens  of  Hertz  to  a  few  hundred  kHz).  In 
order  to  minimize  the  number  of  unknowns  to  invert  for,  the  polarizability  tensor 
is  written  as  a  time-independent  orientation  matrix  multiplied  by  a  time-dependent 
diagonal  intrinsic  polarizability  tensor.  Similarly,  position  is  supposed  to  be  time 
independent  so  that  both  position  and  orientation  angles  are  inverted  only  once 
using  all  time  channels  collected  by  the  instrument.  Moreover,  using  the  dipole 
approximation,  we  are  able  to  compute  the  Jacobian  in  closed  form  for  instruments 
with  either  square  or  circular  primary  field  coils,  thus  contributing  to  the  speed  of 
the  algorithm.  Validating  results  are  shown  based  on  measurement  data  collected 
with  two  EMI  sensors  on  various  types  of  UXO. 

4.2.a  Introduction 

Numerous  military  campaigns  and  warfare  practice  over  the  last  decades  have 
left  millions  of  acres  of  land  contaminated  worldwide  with  unexploded  ordnance 
(UXO),  projectiles  with  explosive  cores  and  metallic  shells  that  did  not  explode 
upon  impact  and  penetrated  into  the  ground.  The  cleaning  of  these  fields  is  cur¬ 
rently  a  daunting  and  hugely  expensive  task  that  consumes  a  budget  of  several 
million  dollars  in  the  United  States  alone.  Such  prohibitively  large  cost  is  due  to 
a  very  large  false-alarm  rate  associated  with  the  detection  of  UXO  with  current 
remote  sensing  technologies.  Differentiating  between  unexploded  and  exploded 
ordnance  is  a  technological  challenge  that  is  yet  to  be  overcome  in  realistic  field 
situations,  where  clutter  and/or  geological  factors  can  significantly  jam  the  sig¬ 
nals  used  by  remote  sensing  algorithms.  As  a  result,  current  cleaning  campaigns 
are  based  on  overly  simple  metal  detection  schemes  with  little  discrimination  in¬ 
volved,  yielding  hugely  time  consuming  and  expensive  campaigns  which  result 
most  of  the  time  in  the  excavation  of  harmless  objects. 

In  parallel,  researchers  have  spent  decades  investigating  technologies  and  re- 
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fining  algorithms  to  detect  and  discriminate  metallic  subsurface  objects  better.  For 
example  radar  sensors  operating  in  the  microwave  regime,  which  have  been  so 
successful  when  pointing  to  the  sky,  have  been  turned  upside  down  to  send  and 
receive  signals  into  and  from  the  ground.  Advanced  signal  processing  techniques 
have  proven  the  usefulness  of  this  approach  [15,  16],  but  also  revealed  its  limita¬ 
tions  due  to  the  variations  in  soil  surfaces  and  volume  compositions  which  quickly 
obscure  the  signatures  of  the  subsurface  targets.  As  an  alternative,  electromagnetic 
induction  (EMI)  which  utilizes  magnetic  fields  at  much  lower  frequencies  (from 
tens  of  Hertz  to  a  few  hundred  kHz)  has  proven  remarkably  robust.  The  associ¬ 
ated  magneto-quasistatic  regime  ensures  that  displacement  currents  are  negligible 
compared  to  conduction  currents  in  the  metallic  targets.  The  soil,  assumed  here 
to  be  non-magnetic  and  non-conductive,  is  transparent  so  that  the  UXO  can  be 
simply  considered  in  free  space.  Various  models  have  been  proposed,  from  the 
simple  but  effective  dipole  model  [17-20],  to  more  elaborated  theoretical  studies 
on  canonical  shapes  [21-23]  as  well  as  more  flexible  numerical  models  [24,  25]. 
The  purpose  of  these  various  techniques  is  to  improve  the  discrimination  capabil¬ 
ities  of  UXO  in  realistic  environments  [26-29]. 

A  major  drawback  of  the  EMI  regime,  however,  is  its  sharp  held  decay  (ap¬ 
proximately  as  1  /r3  at  large  distances  compared  to  the  primary  coil  size,  with  r 
being  the  target- sensor  distance)  which  makes  it  difficult  to  use  EMI  frequencies 
at  ranges  beyond  a  few  times  the  size  of  the  coils  or  when  identifying  multiple 
UXO  in  close  proximity  to  one  another  as  well  as  in  close  proximity  to  large  clut¬ 
ter  items.  For  example,  even  a  small  but  shallow  metal  object  can  overshadow  the 
signal  of  a  deeper  UXO  producing  an  erroneous  classification.  A  related  problem 
occurs  when  targets  of  similar  sizes  are  buried  alongside  one  another  at  compa¬ 
rable  depths  [7].  In  both  situations,  the  detection  and  identification  can  be  con¬ 
siderably  improved  by  considering  a  multi-target  configuration,  whereby  multiple 
time  signatures  are  simultaneously  inverted  for  and  used  in  the  decision  scheme. 
One  of  the  first  attempts  along  this  direction  collaboratively  used  EMI  and  radar 
signals  [9,  10].  This  approach,  however,  revealed  to  be  both  time  consuming  and 
costly  due  to  the  combination  of  two  sets  of  hardware  and  two  analysis  algorithms 
involved.  Other  attempts  concentrated  on  using  only  EMI  data  with  signal  pro¬ 
cessing  techniques  for  sources  separation  [5].  Successful  results  were  reported 
on  synthetic  data  (i.e.  computer  generated)  while  more  validation  is  expected  us¬ 
ing  real  UXO  and  EMI  sensors.  Similarly,  an  independent  component  analysis 
revealed  the  possibility  of  extracting  the  signature  of  a  UXO  in  the  presence  of 
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multiple  clutter  items,  but  the  method  concentrated  on  the  specific  situation  of 
known  UXO,  i.e.  a  known  signature  which  had  to  be  extracted  from  a  sum  of  mul¬ 
tiple  signatures  [30].  More  recently,  positions,  orientations,  and  time  dependent 
polarizabilities  were  inverted  on  two-target  configurations  from  both  synthetic  as 
well  as  sensor  data  [31].  The  method,  based  on  a  multi-start  search  and  separa¬ 
tion  of  the  linear  and  nonlinear  parameters,  was  shown  to  yield  good  agreement 
with  ground  truth.  However,  as  the  authors  point  out,  the  strong  nonlinearity  of 
the  problem  required  a  good  initial  guess,  which  was  provided  by  sampling  a  re¬ 
gion  of  interest.  This  requirement  could  be  alleviated  by  using  a  multiple  signal 
classification  approach  [32,  33]  for  position  estimation  (nonlinear  part)  followed 
by  a  direct  solution  of  the  polarizabilities  (linear  part)  using  either  a  least  square 
method  [34],  a  Kalman  filter  approach  [35],  or  other.  Yet,  this  approach  presents 
some  limitation  when  monostatic  data  are  used  in  the  inversion  process. 

Trying  to  alleviate  some  of  these  pitfalls,  we  propose  in  this  paper  an  algorithm 
significantly  improved  compared  to  our  previous  work  [36]:  factorization  of  the 
orientation  angles  from  the  time  dependent  unknown  vector,  incremental  step  size 
in  the  Gauss-Newton  algorithm,  analytical  computation  of  the  Jacobian  matrix 
for  square  and  circular  primary  held  coils,  and  solution  selection.  The  resulting 
new  algorithm  shows  robustness  to  initial  conditions  that  removes  the  need  of  user 
defined  regions  of  interest,  numerical  stability  that  reduces  the  number  of  realiza¬ 
tions  necessary  to  filter  out  local  minima  solutions,  and  computational  speed.  The 
method  is  here  extensively  validated  using  measured  data  from  two  types  of  EMI 
sensors  (MPV  and  TEMTADS)  on  a  variety  of  UXO. 

4.2.b  Problem  configuration  and  EMI  sensors 

We  consider  the  problem  illustrated  in  Figure  4,  whereby  a  sensor  at  position 
rs  illuminates  multiple  targets  described  by  their  positions  A  =  xxf  +  yyi  +  zzf, 
elevations  Of  and  azimuths  (j)f  where  t  =  1, . . . ,  N  (N  =  3  is  shown  in  the  fig¬ 
ure).  We  suppose  all  N  targets  to  be  bodies  of  revolution,  modeled  as  a  three- 
axis  dipole  with  parameters  ff,  Of  .  (y  and  intrinsic  polarizabilities  (3  f  (see  sub¬ 
sequent  Eq.  (4.5c)).  The  transmitter  and  receiver  sensors  are  incorporated  into 
our  algorithm  to  closely  model  the  hardware  configuration  of  either  the  TEM¬ 
TADS  [37,  38]  or  the  MPV  [39]  instruments,  receiver  sensors  separated  by  40  cm 
(center  to  center)  in  the  (xy)  plane  and  collecting  z  components  of  the  secondary 
magnetic  held  induced  by  the  target.  Each  sensor  is  composed  of  a  square  trans- 
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Figure  4:  Problem  configuration:  an  EMI  sensor  surveys  an  area  where  N  targets  are  buried  (N  —  3 
in  the  figure).  The  unknowns  of  the  problem  are  the  targets  (intrinsic)  polarizabilities  themselves 
as  well  as  their  positions  r#  and  orientations  (0£,</y).  The  soil  is  supposed  to  be  non-conductive 
(<7  =  0)  and  non-permeable  (/ir  =  1). 


mitting  coil  (side  dimension:  35  cm)  and  a  co-centered  square  receiving  coil  (side 
dimension:  25  cm)  which  are  time  decoupled:  during  the  first  half  of  the  cycle, 
the  transmitting  coils  emits  a  primary  field  which  excites  all  metals  in  its  vicin¬ 
ity  whereas  during  the  second  half  of  the  cycle,  the  receiving  coil  measures  the 
derivative  of  the  flux  generated  by  all  ringing  metals.  The  hardware  electronics 
integrates  the  flux  to  yield  a  quantity  proportional  to  the  magnetic  field  which  is 
measured  typically  after  the  first  100  /is  to  let  the  instrument  transients  die  away. 
The  targets  themselves  typically  ring  for  many  milliseconds  so  that  measurements 
are  gathered  until  about  25  ms  logarithmically  sampled  in  120  time  channels. 
Measurements  can  be  performed  in  either  a  monostatic  mode,  in  which  case  all 
the  25  sensors  sequentially  transmit  and  receive  yielding  25  data  points,  or  in  a 
bistatic  mode  in  which  case  each  transmitter  fires  sequentially  and  measurements 
are  gathered  by  all  25  receivers  yielding  625  data  points.  We  model  the  transmit¬ 
ters  as  square  current  loops  whose  primary  field  can  be  written  as  the  addition  of 
four  sections  of  straight  wire  which  all  carry  a  unit  constant  current.  The  primary 
field  from  each  wire  is  given  by 


^wire 


(z\  z\  \  z\ 

r\-r2)-ri2 


(4.1) 
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where  a  hat  denotes  a  unit  vector  and  a  bar  denotes  a  regular  vector,  r\  and  ?2  are 
the  unitary  vectors  between  each  extremities  of  the  straight  wires  and  the  observa¬ 
tion  point,  f\2  =  ri  —  n,  and  s  =  ri2  x  r\.  The  received  signal  at  each  receiver  is 
simply  computed  by  integrating  the  secondary  magnetic  held  over  the  surface  of 
the  corresponding  receiving  square  coil. 

The  MPV  (Man  Portable  Vector)  [39]  sensor  is  very  different  from  the  TEM- 
TADS  in  that  it  is  much  smaller  and  lighter.  It  uses  a  pair  of  coaxial  circular 
transmitting  coils  (of  radius  37.5  cm  and  vertical  separation  15.6  cm)  and  a  series 
of  five  square  receiving  coils  (of  side  dimensions  10  cm)  which  can  measure  all 
components  of  the  secondary  held.  The  primary  held  is  modeled  as  that  produced 
by  a  circular  current  loop  of  radius  a  and  unit  current.  In  the  cylindrical  coordinate 
system,  the  held  at  (p,z)  is  written  as  [18,  40]: 


Hp(p,z) 

Hz(p,z) 


1  1  z 

y/ ( a  +  p)2  +  z 2  P 
1  1 
y/ (a  +  p) 2  +  z2 


a2  +  p2  +  z2 
(p  -  a)2  +  z2 


(p  -a)2+z2 


E{k)  —K{k)  , 
E{h)  +  K{k)  , 


(4.2a) 

(4.2b) 


where  K(k)  and  E(k)  are  the  elliptic  integrals  of  the  hrst  and  second  kind,  re¬ 
spectively,  with  k 2  =  4ap /  [(p  +a)2  +  z2]-  The  received  signal  is  computed  by 
integrating  the  secondary  magnetic  held  over  the  surface  of  the  receiving  perpen¬ 
dicular  square  coils.  Like  for  the  TEMTADS,  measurements  are  gathered  during 
off-time  cycles  (when  the  primary  held  is  turned  off),  between  about  0.4  ms  and 
23  ms  sampled  in  30  time  channels.  Inherent  to  its  construction,  the  operating 
mode  of  the  MPV  is  a  mixture  of  monostatic  and  bistatic  mode,  typically  coined 
multi-static  (the  receivers  are  tied  to  the  transmitter  but  not  collocated).  Advan¬ 
tages  compared  to  the  TEMTADS  are  its  vectorial  property  and  its  maneuverabil¬ 
ity  in  complex  terrain  whereas  its  disadvantages  are  its  smaller  held  of  view  and 
the  necessity  of  tracking  its  position  and  orientation  during  data  acquisition. 


4.2.c  Data  analysis 

The  fast  held  decay  in  the  EMI  regime  ensures  that  near  held  effects  are  conhned 
to  a  very  small  volume  surrounding  the  buried  target.  As  a  consequence,  the  sec¬ 
ondary  held  at  the  sensor  location  is  almost  identical  to  that  of  a  triaxial  dipole, 
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for  which  the  magnetic  held  H  is  given  by  [18]: 


H(rs) 


3  rts{rts-m) 


(4.3) 


where  m  =  xmx  +ymy  +z,mz  is  the  dipole  moment  of  the  target,  rt  is  the  target  loca¬ 
tion,  rs  is  the  sensor  location,  rts  =  rs  —  rt  and  rts  =  |  \rts\  |.  It  should  be  mentioned 
that  the  dipole  analysis  remains  an  approximation  which  has  proven  accurate  in 
the  overwhelming  majority  of  UXO  configurations.  However,  its  limitation  is  par¬ 
ticularly  obvious  with  heterogeneous  targets  or  when  multiple  targets  are  in  very 
close  proximity  to  one  another,  typically  within  a  sub-centimeter  distance,  with 
one  or  both  targets  being  ferrous  [41].  In  this  situation,  coupling  between  eddy 
currents  of  the  two  objects  can  happen  at  a  level  sufficient  to  alter  the  secondary 
magnetic  held  at  the  receiver  location  [42].  Such  conhguration  is,  however,  more 
theoretical  than  practical  and  to  our  knowledge,  has  not  been  encountered  beyond 
laboratory  tests.  As  a  consequence,  we  suppose  here  that  the  separation  between 
the  targets  is  sufficient  to  ensure  that  the  various  UXO  respond  as  independent 
dipoles  with  no  coupling  effects  [43].  Within  such  independent  scattering  approx¬ 
imation,  the  magnetic  held  of  N  targets  located  at  ?i  (f  =  1 , . . . ,  N)  and  with  dipole 
moments  mi  is  simply  obtained  as  the  algebraic  sum  of  Eq.  (4.3): 


H(rs) 


N 


I 


1 

4*4 


3 ris[rh  ■  mi) 


(4.4) 


Note  that  we  assume  N  to  be  known  in  Eq.  (4.4),  which  is  an  assumption  further 
discussed  in  Section  4.2.e.  The  dipole  moment  fhi  is  obtained  from  the  self¬ 
polarizability  Pi  of  the  target  along  three  directions,  weighted  by  the  primary  held 
at  the  target  location.  If,  in  addition,  the  target  is  rotated  with  respect  to  the  global 
reference  frame,  the  rotation  angles  need  to  be  incorporated.  Restricting  ourselves 
to  bodies  of  revolution  as  most  UXO  are,  the  dipole  moment  mi  can  be  written  as 

(4.5a) 


where 


(cos  0i  cos  ( pi  cos  0i  sin  (j)£ 
—  sin  (j>i  cos  <j)i 

sin  6i  cos  (j>i  sin  6i  sin  (j)i 


(4.5b) 
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is  the  rotation  matrix  with  the  angles  Of  and  (j)f, 

=  (K  o  o\ 

Pe  =  o  Pyt  0  ,  (4.5c) 

\  0  0  pj 

where  ((3Xn Pyn  (5Ze)  are  the  time-dependent  polarizabilities  along  the  three  direc¬ 
tions  of  the  reference  frame  local  to  the  target,  and  Hpr{?f)  is  the  primary  field 

at  the  target  location.  Note  that  two  components  of  /3  ^  are  expected  to  be  iden¬ 
tical  due  to  the  BOR  assumption,  even  though  this  property  is  not  enforced  in 
Eq.  (4.5c).  For  example,  in  a  bistatic  configuration,  the  transmitter  remains  at  a 
unique  location  while  the  receiver  spans  a  series  of  measurement  points.  As  a 
consequence,  the  primary  field  Hpr{?f)  is  unique  and  so  is  fhf.  In  a  monostatic 
configuration,  however,  where  the  transmitter  and  receivers  are  collocated  and 
span  together  the  measurement  points,  Hpr{?f)  varies  from  point  to  point,  and  so 
does  fhf.  The  polarizability  tensor  (3  p,  however,  is  invariant  and  can  be  used  by  a 
classifier  for  each  target. 

Our  inversion  algorithm  is  based  on  a  Gauss-Newton  solution  of  the  normal 
equation  assuming  N  targets  present  in  the  field  of  view  of  the  sensor.  The  simulta¬ 
neous  inversion  of  the  time-independent  position,  orientation,  and  time-dependent 
polarizabilities  of  the  targets  is  expressed  as  a  least-squares  minimization  between 
the  measured  field  and  the  predicted  field  computed  from  Eq.  (4.4).  Hence,  the 
unknown  vector  is  written  as 


x  =[ai . . .  . . .  aN  bi . . . be . . .  bN]T , 

*£  =[xtytzi  Of  (j)f\T , 


(4.6a) 

(4.6b) 

(4.6c) 


where  tj  (/  =  1  denotes  the  time  channels  and  where  the  subscript  ‘T’ 

denotes  the  transpose.  Eq.  (4.6c)  indicates  that  the  algorithm  solves  for  a  single 
position  and  a  single  orientation  per  target,  with  multiple  polarizabilities  across 
the  time  history  of  the  transmitter  for  all  targets  simultaneously.  Note  that  it  is 
also  possible  to  invert  for  a  more  general  polarizability  tensor  with  six  indepen¬ 
dent  elements  [31,  44,  45].  Although  such  approach  has  been  shown  to  be  valid 
and  accurate,  it  is  not  adopted  here:  inverting  a  full  tensor  at  each  time  channel  is 
analogous  to  inverting  for  the  rotation  angles  at  each  time  channel.  Within  a  good 
approximation,  however,  these  angles  can  be  considered  time  independent  so  that 
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we  factor  them  out  of  the  polarizability  tensor  and  included  them  in  the  time  in¬ 
dependent  part  of  the  minimization  vector,  thus  drastically  reducing  the  number 
of  unknowns  to  invert  for.  Note  that  similarly,  the  algorithm  inverts  for  a  unique 
position  rather  than  a  time-dependent  position  at  every  time  channel.  Consider¬ 
ing  the  position  and  angles  to  be  time  independent  is  an  approximation  justified 
in  practice.  Indeed,  although  various  parts  of  the  UXO  may  respond  differently 
in  time,  effectively  yielding  slightly  different  positions  and  orientations  of  the  tri- 
axial  dipole  used  to  best  match  the  measured  secondary  held,  we  have  found  this 
effect  to  be  small  in  the  cases  considered  in  this  paper:  inverting  the  positions  at 
every  time  channels  yields  similar  estimates  except  at  later  times  when  the  held 
values  reach  the  noise  level  (data  not  shown). 

The  task  of  the  minimization  algorithm  is  to  estimate  the  unknown  parameter 
x  so  as  to  minimize  the  sum  of  squares 

minis'  =  min  (A Hi  AHS)  (4.7) 

where  A Hz  =  (/7mca  —  Hz)  with  /7mca  and  Hz  being  the  measured  and  computed 
z  components  of  the  magnetic  held,  respectively.  If  needed,  the  addition  of  other 
components  of  the  magnetic  held  to  the  minimization  problem  is  straightforward 
(e.g.  the  TEMTADS  measures  only  Hz  whereas  the  MPV  measures  all  three 
components).  The  corresponding  normal  equation  is  written  as 

7zT7z  Ax  =  7zT  AHZ  (4.8) 

where  Ax  =  x'+1  —  x*  is  the  unknown  vector  which  provides  the  updated  quantities 
of  x  at  each  iteration,  i  is  the  iteration  number,  Jz  is  the  Jacobian  matrix.  The 
Gauss-Newton  update  equation  is  therefore 

x/+1  =  x;'  +  Ax  =  x*  +  Qz  7Z)  ~ 1 7Z  ■  (/7f ea  -  Hz)  (4.9) 

The  vector  notation  used  on  /7mca  and  Hz  refers  to  the  fact  that  the  magnetic 
held  is  measured  and  computed  over  a  series  of  K  observation  points.  The  Jaco¬ 
bian  matrix  is  therefore  a  block  matrix  of  size  [Ntc  x  K.N(5  +  3 Ntc)]  where  Ntc  is 
the  number  of  time  channels.  Splitting  the  Jacobian  matrix  into  time  independent 
and  time  dependent  parts,  we  write: 


(4.10) 
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where  Az,  of  size  [Ntc  x  K,5N],  gathers  the  derivatives  of  the  computed  Hz  held 
with  respect  to  the  three  positions  and  two  angles  across  all  observation  points 
and  time  channels.  Assuming  N  =  1  for  the  sake  of  simplicity  of  notation  in  this 
example: 


(t) 

/ Atj  \ 

4  = 

A’z 

with  A t  = 

^ ze 

(4.11a) 

=* Ntc 

\Az  / 

..  c-A 

and 


A]  — 

^ ze  — 


d  d  d  d  d 


Hhn) 


(4.11b) 


dx£  dyi  d  Z£  d6{  d  (j)L 

where  Hz  {?(,■)  is  the  z  component  of  the  magnetic  held  obtained  from  Eqs.  (4.4) 
and  (4.5)  at  time  channel  tj.  The  matrix  B.  is  a  block  diagonal  matrix  where  each 
block  has  a  size  [Ntc  x  K ,  3Ntc]  (still  under  the  assumption  of  N  =  1)  and  gathers 
the  derivatives  of  the  computed  Hz  held  with  respect  to  the  three  polarizabilities: 


Bz  =  diag[5 ' 


=ti 


z  ’ 


...,5 


z  ’ 


.5 


Ntc 


with  B  -  = 


( Btj  \ 

Dz  1 


jJi 

BZ( 


\4/ 


(4.12a) 


and 


B%  = 


iri) 


(4.12b) 


dfixi  d(5yt  dpze  j 

The  formulation  outlined  above  is  straightforward  to  generalize  to  the  case  of 
N  >  1  as  well  as  when  all  components  of  the  magnetic  held  are  known  (for  ex¬ 
ample  if  the  MPV  sensor  is  used)  by  simply  stacking  the  Jacobian  matrices  and 
right-hand  side  vectors  corresponding  to  other  components,  and  generalizing  the 
unknown  vector  x.  Note  that  all  the  derivatives  are  straightforward  to  compute  ex¬ 
cept  those  with  respect  to  position  since  they  require  the  derivative  of  the  primary 
held  as  well,  in  view  of  Eq.  (4.5a).  This  task  can  still  be  performed  analytically 
for  both  the  TEMTADS  and  MPV  sensors  using  Eqs.  (4.1)  and  (4.2)  and  a  pre- 
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cise  knowledge  of  the  sensors’  geometries,  which  contributes  to  the  computational 
efficiency  of  the  algorithm. 

The  actual  implementation  of  Eq.  (4.9)  is  slightly  modified  by  the  addition  of  a 
weighting  factor  a  <  1  in  order  to  control  the  step  size: 


xi+1  =x‘  +  a(JzT  Jz)  1 JZT  ■  (H™ :a  -  H.) 


mea 


(4.13) 


Empirical  investigations  revealed  an  increasing  importance  of  a  as  the  number 
of  targets  to  invert  for  increased:  from  practically  unnecessary  for  N  =  1  to  im¬ 
portant  for  N  >  3.  While  further  analytical  investigations  might  be  necessary  to 
determine  the  optimal  a  for  each  configurations,  results  reported  hereafter  have 
been  obtained  with  a  =  0.5,  which  appears  to  offer  a  good  compromise  between 
accuracy  and  convergence  speed. 

4.2. d  Inversion  results 

Despite  an  extensive  validation  campaign  on  synthetic  data,  we  concentrate  and 
only  present  here  results  from  data  that  have  been  acquired  by  real  sensors  on 
real  UXO.  Data  have  been  collected  using  both  the  TEMTADS  [46]  and  MPV 
instruments  [47]  and  the  inversion  results  are  presented  in  the  following  two  sec¬ 
tions  on  various  multi-target  configurations.  It  is  important  to  keep  in  mind  that 
measurements  are  performed  on  real  UXO  which  therefore  have  certain  physical 
dimensions,  whereas  the  algorithm  is  inverting  for  one  tri-axial  dipole  per  UXO. 
The  location  of  these  equivalent  dipoles  is  not  precisely  defined:  it  is  close  to  the 
center  of  the  target  if  the  latter  is  buried  horizontally  (i.e.  parallel  to  the  sensor), 
whereas  it  migrates  toward  the  tip  for  vertical  targets  where  one  extremity  is  closer 
to  the  sensor  than  the  other.  The  positions  also  depend  on  the  composition,  geom¬ 
etry,  and  external  features  of  each  target  such  as  fins,  conductive  rings,  etc.  The 
exact  location  at  which  a  single  tri-axial  dipole  should  be  located  in  order  to  best 
reproduce  the  secondary  held  of  a  given  UXO  is  therefore  not  expected  to  match 
exactly  the  positions  recorded  at  the  data  acquisition  stage  (which  correspond  typ¬ 
ically  to  the  exact  center  or  exact  tip  of  the  UXO).  The  accuracy  of  the  subsequent 
results  should  be  judged  with  these  considerations  in  mind. 

4.2. d.(l)  MPV  data  The  MPV  is  a  vectorial  multi- static  EMI  instrument  which 
collects  data  as  it  is  waved  above  the  surface.  Due  to  its  intrinsic  portable  nature, 
the  trajectory  of  the  MPV  as  well  as  its  orientation  (e.g.  pitch,  roll,  and  yaw)  need 
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to  be  tracked  during  the  entire  data  acquisition  time.  The  data  set  thus  collected 
is  then  converted  into  components  of  the  magnetic  held  along  three  pre-dehned 
directions,  and  subsequently  processed  by  our  algorithm.  The  results  reported 
hereafter  have  been  obtained  using  the  vectorial  information  of  all  five  sensors, 
even  though  more  in  depth  studies  showed  that  not  all  carry  a  critical  information 
necessary  for  proper  inversion  (data  not  shown). 

A  first  series  of  validating  tests  of  our  MPV  model  was  performed  on  single 
targets  and  yielded  good  results  in  all  the  cases  considered  for  both  the  positions, 
orientations,  and  time  dependent  polarizabilities.  Table  3  lists  a  few  cases  with 
the  corresponding  inverted  positions  and  orientations  which  show  a  reasonable 
agreement  with  the  ground  truth.  We  recall  here  that  the  ground  truth  values  typ¬ 
ically  correspond  to  distances  to  the  center  of  the  UXO,  which  might  be  a  few 
centimeters  deeper  from  the  optimized  location  of  an  equivalent  dipole  (nearer  to 
the  sensor).  For  example,  the  elevation  of  the  60-mm  UXO  was  set  to  0  =  270 
degrees,  which  corresponds  to  the  target  being  horizontal.  The  exact  location  of 
the  principal  responding  magnetic  dipole  of  the  target  may  not  correspond  to  the 
geometrical  center  of  the  actual  target.  Therefore,  a  first  order  correction  to  the 
inverted  depth  for  this  case  to  account  for  the  radius  of  the  target,  i.e.  30  mm, 
augments  the  inverted  depth  from  48  cm  to  51  cm,  in  good  agreement  with  the 
ground  truth.  Corrections  to  other  targets  are  not  as  straightforward  to  incorporate 
because  of  the  non-trivial  elevation  angles.  The  largest  deviation  in  depth  amounts 
to  10  cm  and  occurs  with  the  105-mm  with  0  =  — 19  degrees.  The  small  elevation 
indicates  that  the  target  is  close  to  vertical  so  that  an  important  deviation  between 
the  center  position  and  the  location  of  the  equivalent  dipole  is  expected,  especially 
because  this  target  is  very  long.  In  addition  to  the  mismatch  between  the  center 
of  the  target  and  the  location  of  its  equivalent  dipole,  part  of  the  discrepancies 
in  Table  3  can  also  be  attributed  to  measurement  and  instrument  noise  while  yet 
another  part  could  be  due  to  algorithmic  inaccuracies  (for  example  due  to  a  too 
small  data  set  to  optimize  the  unknown  vector  to  a  desired  level  of  accuracy).  The 
overall  good  results,  however,  suggest  that  the  algorithm  is  effective  at  inverting 
MPV  data.  This  conclusion  is  also  supported  by  the  time  dependent  polarizabili¬ 
ties  obtained  for  these  cases,  a  subset  of  which  is  shown  in  Figure  5.  As  expected, 
the  81 -mm  UXO  shows  one  polarizability  element  consistently  stronger  than  the 
other  two  which  are  almost  identical  and  reveal  a  geometric  azimuthal  symmetry. 
The  BLU26  exhibits  similar  characteristics  but  the  three  polarizabilities  are  closer 
to  each  other  due  to  the  close  to  spherical  shape  of  this  UXO  (see  inset). 
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Target 

Length 

[mm] 

X 

Ground  truth 
[cm]  [degrees] 

y  z  0  (j) 

X 

Inverted  results 
[cm]  [degrees] 

y  z  0  (j) 

A1  sphere 

0 

22 

26 

- 

- 

0 

22 

29 

- 

- 

Brass  sphere 

0 

22 

26 

- 

- 

0 

22 

28 

- 

- 

Steel  sphere 

0 

22 

26 

- 

- 

0 

22 

28 

- 

- 

81 -mm 

490 

-23 

22 

52 

-18 

0 

-14 

22 

50 

-20 

-10 

105 -mm 

642 

-20 

22 

63 

-19 

180 

-12 

23 

50 

-30 

176 

BLU26 

70.31 

0 

22 

38 

0 

0 

0 

22 

34 

- 

- 

57-mm 

170.47 

5 

22 

55 

307 

180 

4 

22 

51 

310 

160 

60-mm 

245 

0 

22 

51 

270 

0 

0 

22 

48 

260 

7 

Table  3:  Inverted  positions  of  single  targets  using  data  from  the  MPV  instrument.  The  time  sig¬ 
natures  for  the  81 -mm  and  the  BLU26  are  shown  in  Figure  5.  Positions  are  rounded  to  the  closest 
integer. 


Time  [ms] 


Figure  5:  Time  domain  polarizabilities  for  two  types  of  UXO  (81 -mm  and  BLU26)  measured  and 
inverted  separately  from  MPV  data. 
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Two-target  results  present  mitigated  yet  consistent  results  within  the  data  set 
available  to  us.Two  configurations  were  analyzed:  that  of  two  40-mm  UXO  at 
the  same  depth,  and  that  of  a  40-mm  along  with  81 -mm  UXO  at  varying  depths, 
both  cases  for  different  lateral  separations.  Inverted  positions,  summarized  in  Ta¬ 
ble  4,  indicated  a  good  performance  when  the  targets  are  at  the  same  depth,  with 
yet  improved  accuracy  for  larger  lateral  separations  as  intuitively  expected.  This 
is  evident  from  the  inverted  polarizabilities  of  the  two  40-mm  targets,  shown  in 
Figure  6(a):  while  the  four  strong  components  (two  targets  in  two  configurations) 
are  almost  identical  in  both  cases,  one  of  the  weak  components  is  slightly  under¬ 
estimated  for  dx= 25  cm  (circle  marks).  The  good  performance  at  identical  depth 
is  somewhat  maintained  with  the  two  different  targets,  despite  a  small  error  in  the 
inverted  lateral  separation  which  can  be  mostly  attributed  to  the  physical  width  of 
the  targets.  Incorporating  this  correction  of  about  6  cm,  the  inverted  separation 
would  be  about  21  cm  instead  of  the  measured  25  cm.  The  overall  results  worsen, 
however,  when  the  two  targets  are  at  different  depths,  with  zi  in  Table  4  often 
under-estimated.  The  polarizabilities,  shown  in  Figure  6(b)  reflect  this  behavior: 
while  the  larger  target  is  consistently  well  predicted,  the  weaker  one  is  better  pre¬ 
dicted  at  similar  depth  (see  the  dipping  weak  component  for  dz= 25  cm,  indicating 
a  sign  change  and  thus  a  weak  corresponding  component  of  the  secondary  mag¬ 
netic  held). 

In  depth  investigations  of  these  cases  suggested  that  the  data  set  itself  might  not 
be  totally  consistent  from  grid  point  to  grid  point,  but  contains  variations  that  are 
likely  at  the  origin  of  some  of  the  errors  in  the  inversion  results.  As  a  matter  of  fact, 
although  the  natural  operating  mode  of  the  MPV  is  to  be  waved  above  the  targets, 
data  provided  to  us  has  been  collected  in  an  inverted  configuration,  whereby  the 
MPV  was  held  fixed  and  the  targets  were  moved  on  a  5  x  5  square  grid  with  20  cm 
separation  between  adjacent  points.  Repeatability  between  grid  point  therefore 
cannot  be  guaranteed  since  the  two  targets  can  easily  move  with  respect  to  each 
other,  either  in  distance  or  in  orientation.  Since  our  algorithm  assumes  an  invariant 
configuration  from  grid  point  to  grid  point  (akin  to  when  the  targets  are  buried), 
these  errors  would  induce  inversion  inaccuracies.  Yet,  the  overall  good  agreement 
between  the  inverted  results  and  reported  ground  truth  confirm  the  possibility  of 
using  the  MPV  for  the  detection  and  classification  of  multi-targets. 


4.2.d.(2)  TEMTADS  data  measurement  campaign  using  the  TEMTADS  sensor 
has  been  undertaken  by  the  Naval  Research  Laboratory  on  one,  two,  as  well  as 
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Time  [s] 


(a)  Two  40-mm  UXO  inverted  in  two  configurations:  (dx,dz)=( 25,0)  and 
(40,0)  cm.  The  visibly  weaker  components  (circle  marks)  corresponds  to 
dx= 25  cm. 


Time  [s] 

(b)  81-mm  and  40-mm  UXO  inverted  in  two  conhgurations: 

(dx\dz)=(25S))  and  (25,25)  cm.  The  weaker  and  dipping  component  (di¬ 
amond  marks)  corresponds  to  dz= 25  cm. 

Figure  6:  Time  domain  polarizabilities  for  two  two-target  configurations  measured  by  the  MPV. 
Inverted  positions  are  summarized  in  Table  4. 
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UXO 

X\ 

yi 

Truth  [cm] 

zi  x2 

>’2 

Z2 

Xl 

yi 

Inverted  [cm] 

Zi  x2 

J2 

Z2 

40-mm  &  40-mm 

0 

0 

40 

25 

0 

40 

2 

-l 

40 

23 

0 

38 

40-mm  &  40-mm 

0 

0 

40 

40 

0 

40 

1 

-l 

39 

40 

-2 

39 

8 1  -mm  &  40-mm 

0 

0 

40 

25 

0 

40 

0 

0 

42 

18 

0 

40 

81-mm  &  40-mm 

0 

0 

65 

25 

0 

40 

0 

0 

43 

22 

0 

38 

Table  4:  Inversion  results  on  two-target  configurations  from  data  collected  by  the  MPV  sensor. 
Positions  are  rounded  to  the  closest  integer. 


three-target  configurations.  The  targets  of  interest  were  a  4.2-inch  mortar,  a  nose 
cone,  and  a  baseplate  whose  time  dependent  polarizabilities  are  represented  in 
Figure  7  when  they  are  inverted  as  separate  targets.  Like  the  81-mm  previously, 
the  mortar  exhibits  a  typical  signature  of  a  UXO,  whereby  one  polarizability  dom¬ 
inates  the  other  two  which  are  almost  identical,  indicating  that  the  target  is  a  body 
of  revolution  with  a  preferential  direction.  The  nose  cone  produces  a  signal  that 
is  about  two  orders  of  magnitude  inferior  to  the  mortar,  with  a  mostly  log-log 
linear  time  decay  indicating  a  strong  ferrous  component.  In  addition,  the  three  po¬ 
larizabilities  are  comparable  which  suggests  a  somewhat  homogeneous  shape  in 
three  directions,  akin  to  a  sphere.  The  baseplate  produces  a  signal  of  comparable 
magnitude  to  the  mortar  with  the  strongest  dipole  switching  in  time  from  the  x 
direction  to  the  z  direction.  These  signatures  have  been  obtained  using  the  method 
presented  in  the  previous  section  and  have  been  independently  validated  [48,  49]. 
The  inverted  positions,  listed  in  the  captions  of  Figure  7,  match  the  true  values 
with  a  centimeter  accuracy  and  validate  our  algorithm  as  well  as  the  associated 
TEMTADS  model  in  the  case  of  single  targets.  We  emphasize  also  the  fact  that 
the  inverted  angles  are  typically  within  a  few  degrees  of  their  expected  values. 

Two-target  configurations  consisted  of  combining  the  mortar  with  one  of  the 
other  two  targets  and  varying  their  respective  lateral  distances  while  keeping  their 
depths  constant.  The  lateral  distances  were  varied  from  x  =  0  cm  to  x  =  100  cm  at 
10  cm  intervals  (yielding  1 1  configurations).  At  x  =  0  the  targets  are  immediately 
on  top  of  each  other,  which  is  known  to  be  a  challenging  configuration  to  invert 
due  to  the  strong  EMI  signal  decay.  On  the  other  hand,  the  position  x  =  0  is 
advantageous  because  the  target  is  surrounded  by  sensors  and  therefore  data  carry 
information  from  multiple  views  around  the  targets.  At  x  =  100  cm,  however, 
the  target  is  at  the  edge  of  the  instrument  and  not  immediately  under  a  receiver. 
The  amount  of  quality  data  is  therefore  limited  in  this  case:  multi-looks  around 
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Time  [ms] 


(a)  Mortar.  Inverted  positions:  (-0.5, 0.6, -60.8)  [cm].  True  positions:  (0,0,60)  [cm]. 


(b)  Nose  cone.  Inverted  positions:  (0.2,0. 2, -30. 2)  [cm].  True  positions:  (0,0,30)  [cm]. 


(c)  Baseplate.  Inverted  positions:  (0.4,0. 1, -47.6)  [cm].  True  positions:  (0,0,49)  [cm]. 


Figure  7:  Time  domain  polarizabilities  and  inverted  positions  for  three  independent  UXO  using 
TEMTADS  data. 
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the  target  are  not  available  and  many  receivers  are  too  far  away  to  retain  usable 
information.  As  a  result,  inversion  results  are  expected  to  be  best  away  from  the 
two  extreme  positions  x  =  0  and  x  =  100  cm. 

A  systematic  inversion  of  all  1 1  configurations  was  carried  out  on  both  sets  of 
measurements  (mortar  +  nose  cone  and  mortar  +  baseplate).  The  inverted  po¬ 
sitions  for  the  first  set  are  listed  in  Table  5  and  the  time  signatures  are  shown  in 
Figure  8(a).  The  results  show  polarizabilities  of  the  mortar  practically  identical  for 
all  1 1  configurations,  as  expected  since  the  polarizabilities  (3  are  an  intrinsic  prop¬ 
erty  of  the  target,  independent  of  its  location  and  orientation.  The  polarizabilities 
of  the  nose  cone  present  more  variations,  primarily  due  to  the  fact  that  some  con¬ 
figurations  correspond  to  a  nose  cone  located  close  to  either  x  =  0  and  x  =  100  cm. 
Nonetheless,  all  the  signatures  are  seen  to  be  linear  and  in  good  relationship  with 
one  another.  The  inverted  positions  of  both  targets  remain  within  ±  1  cm  accuracy 
for  x  between  20  cm  and  60  cm,  except  for  the  depth  of  the  nose  cone  whose  error 
can  amount  to  a  few  centimeters.  The  largest  error,  of  about  9  cm,  corresponds  to 
the  case  when  the  lateral  offset  is  x  =  100  cm  which  brings  the  target  at  the  edge 
of  the  instrument. 

The  results  for  the  second  set  of  measurements  (mortar  +  baseplate),  summa¬ 
rized  in  Table  6  and  Figure  8(b),  present  more  variations  between  the  1 1  config¬ 
urations  (note  that  the  polarizabilities  of  the  baseplate  have  been  scaled  by  10-2 
in  order  to  separate  them  from  those  of  the  mortar).  The  signature  of  the  mortar, 
although  still  clearly  separated  into  a  strong  component  and  two  smaller  identi¬ 
cal  ones,  exhibits  more  fluctuations  than  in  the  previous  case  but  not  substantial 
enough  to  hinder  a  proper  identification  of  this  target.  The  polarizabilities  of  the 
baseplate  also  show  a  flipping  dipole  behavior  between  x  in  early  time  and  z  in 
late  time,  albeit  somewhat  smoothed  out  in  the  fluctuations  between  all  the  11 
configurations.  Note  again  that  fluctuations  in  the  signature  of  the  second  object 
are  expected  due  to  its  varying  position  which  reaches  the  extremity  of  the  sensor 
in  some  cases. 


A  series  of  eleven  three-target  configurations  were  analyzed  whereby  two  tar¬ 
gets  were  fixed  and  one  was  moving  from  case  to  case.  The  stationary  targets  were 
the  mortar  held  horizontally  at  (x,  y,  z)  =  (0, 0, 60)  [cm]  and  the  baseplate  held  hor¬ 
izontally  at  (x,y,z)  =  (50,0,49)  [cm].  The  moving  target  was  the  nose  cone  held 
vertically  at  (x,y,z)  =  (Jx,0, 30)  [cm]  with  dx  varying  from  0  to  100  cm  in  incre¬ 
ments  of  10  cm.  The  recovered  positions  for  these  configurations  are  summarized 
in  Table  7  and  show  an  overall  good  agreement  with  the  ground  truth.  The  two 
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(a)  Mortar  +  Nose  cone. 
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(b)  Mortal'  +  Baseplate. 


Figure  8:  Time  dependent  polarizabilities  for  two-target  inversions  from  TEMTADS  data.  The 
various  cases  correspond  to  various  lateral  displacements  along  x,  of  the  shallower  target  while  the 
deeper  one  (the  mortar)  is  kept  at  a  constant  location.  True  and  inverted  positions  are  given  in 
Table  5  and  Table  6.  Note  that  the  polarizabilities  of  the  baseplate  in  (b)  have  been  scaled  by  10  2 
in  order  to  visually  separate  them  from  those  of  the  mortar. 
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Mortar 

x  y  z 

Baseplate 

x  y  z 

Truth 

0  0  60 

dx  0  30 

Case 

1 

-1  1  62 

o  2  0  33 

2 

0  1  60 

10  15  0  30 

3 

1  1  61 

20  2  1  0  34 

4 

1  1  61 

30  31  0  34 

5 

1  1  60 

40  41  -1  34 

6 

0  1  61 

50  51  0  32 

7 

1  1  60 

60  6  0  0  33 

8 

0  1  60 

70  67  0  31 

9 

0  1  60 

80  75  0  31 

10 

0  1  60 

90  7  9  0  33 

11 

0  1  60 

loo  80  0  39 

Table  5:  True  and  inverted  positions  of  two  simultaneous  targets  whose  time  dependent  polariz¬ 
abilities  are  shown  in  Figure  8(a)  (TEMTADS  data).  Ground  truth  corresponds  to  a  fixed  mortar 
and  a  moving  nose  cone  along  the  x  axis  by  the  indicated  dx  while  keeping  y  and  z  constants. 
Deviations  of  the  first  few  cases  are  due  to  the  difficulty  in  inverting  accurate  positions  when  the 
two  targets  are  on  top  of  each  other  whereas  deviations  of  the  last  few  cases  are  explained  by  the 
fact  that  the  target  location  extends  beyond  the  physical  size  of  the  sensor.  All  dimensions  are  in 
centimeters  and  positions  are  rounded  to  the  closest  integer. 


stationary  targets  are  well  located  in  the  (xy)  plane  and  the  depths  are  typically 
estimated  within  a  centimeter  for  the  mortar  and  within  a  couple  of  centimeters 
for  the  baseplate,  except  again  in  a  few  isolated  cases  that  correspond  for  the  most 
part  to  targets  at  the  edge  of  the  instrument.  The  positions  of  the  nose  cone  present 
more  fluctuations  while  still  following  the  ground  truth  reasonably  well  in  the  ma¬ 
jority  of  cases.  The  most  problematic  positions  are  still  those  corresponding  to  the 
target  being  at  the  limit  of  the  held  of  view  of  the  instrument  where  an  error  up  to 
9  cm  in  lateral  position  and  8  cm  in  depth  is  recorded.  The  time-dependent  polar¬ 
izabilities,  shown  in  Figure  9,  confirm  these  conclusions.  The  mortar  as  well  as  the 
baseplate  present  very  consistent  signatures  from  case  to  case.  Interestingly,  the 
signatures  of  the  baseplate  are  more  consistent  than  in  the  two-target  inversions  of 
Figure  8(b),  primarily  due  to  the  larger  lateral  separation  with  the  mortar  (50  cm 
in  this  case).  The  polarizabilities  of  the  nose  cone  are  a  few  orders  of  magnitude 
lower  than  those  of  the  mortar  and  baseplate,  but  with  much  more  variations  from 
case  to  case.  We  again  attribute  these  to  two  main  reasons:  overlapping  targets 
and  edge  effects.  Edge  effects  refer  to  the  effect  already  witnessed,  by  which  the 
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Mortar 

x  y  z 

Baseplate 

x  y 

z 

Truth 

0 

0 

60 

dx 

0 

49 

Case 

1 

-4 

1 

53 

0 

-14 

3 

44 

2 

10 

1 

58 

10 

11 

1 

48 

3 

-1 

1 

56 

20 

26 

1 

47 

4 

1 

1 

58 

30 

32 

1 

48 

5 

1 

0 

59 

40 

41 

1 

48 

6 

1 

0 

59 

50 

50 

0 

48 

7 

1 

0 

60 

60 

59 

0 

47 

8 

1 

0 

60 

70 

67 

0 

47 

9 

1 

0 

60 

80 

73 

0 

47 

10 

0 

0 

60 

90 

80 

0 

49 

11 

0 

0 

60 

100 

89 

0 

51 

Table  6:  True  and  inverted  positions  of  two  simultaneous  targets  whose  time  dependent  polar¬ 
izabilities  are  shown  in  Figure  8(b)  (TEMTADS  data).  The  ground  truth  corresponds  to  a  fixed 
mortar  and  a  moving  baseplate  along  the  x  axis  by  the  indicated  dx  while  keeping  y  and  z  constants. 
Deviations  of  the  first  few  cases  are  due  to  the  difficulty  in  inverting  accurate  positions  when  the 
two  targets  are  on  top  of  each  other  whereas  deviations  of  the  last  few  cases  are  explained  by  the 
fact  that  the  target  location  extends  beyond  the  physical  size  of  the  sensor.  All  dimensions  are  in 
centimeters  and  positions  are  rounded  to  the  closest  integer. 


inverted  position  and  time  domain  polarizabilities  of  a  target  located  at  the  edge 
of  the  instrument  can  present  large  errors.  It  is  therefore  expected  that  the  region 
close  to  x  =  100  cm  is  a  source  of  inaccurate  polarizabilities.  The  overlapping 
target  effect  refers  to  the  fact  that  EMI  signatures  are  intrinsically  challenging  to 
separate  for  targets  in  proximity  and  on  top  of  each  other.  In  this  particular  set  of 
measurements,  the  nose  cone  travels  from  x  =  0  to  x  =  100  cm  and  is  directly  on 
top  of  targets  at  x  =  0  (the  mortar)  and  x  =  50  cm  (the  baseplate).  A  case  by  case 
examination  reveals  that  the  position  x  =  0  is  well  inverted  and  that  the  polariz¬ 
abilities  exhibit  the  expected  linear  decay  (on  a  log-log  scale)  visible  for  example 
in  Figure  8(a).  The  polarizabilities  when  the  nose-cone  is  close  to  x  =  50  cm, 
however,  present  much  more  variation  and  can  be  seen  in  Figure  9(a).  The  major 
difference  between  the  two  situations  are  the  distances  between  the  targets  when 
they  coincide  laterally:  30  cm  in  one  case  and  about  20  cm  in  the  other.  It  is 
expected  that  diversification  of  the  positions  by  allowing  for  variations  in  the  y 
direction  would  significantly  improve  the  results.  Note  also  that  in  all  cases,  in¬ 
cluding  those  challenging  ones,  the  inverted  positions  remain  within  reasonable 
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0 

43 
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- 

- 

- 

Table  7:  True  and  inverted  positions  of  three  simultaneous  targets  whose  time  dependent  polariz¬ 
abilities  are  shown  in  Figure  9  (TEMTADS  data).  The  ground  truth  corresponds  to  a  fixed  mortar, 
a  fixed  baseplate,  and  a  moving  nosecone  along  the  x  axis  by  the  indicated  dx  while  keeping  y  and 
z  constants.  All  dimensions  are  in  centimeters  and  positions  are  rounded  to  the  closest  integer. 


bounds. 

4.2.e  Discussion  and  Conclusions 

The  Gauss-Newton  based  algorithm  presented  in  this  paper  has  been  shown  to  be 
effective  at  simultaneously  inverting  for  the  positions,  orientations,  and  time  de¬ 
pendent  polarizabilities  of  multiple  UXO  within  the  held  of  view  of  EMI  sensors. 
Part  of  the  efficiency  is  due  to  the  construction  of  the  unknown  vector,  in  which 
positions  and  orientations  of  the  targets’  equivalent  dipoles  are  assumed  time  in¬ 
dependent  and  inverted  for  only  once  using  the  information  of  all  time  channels. 
Another  important  aspect  of  the  algorithm  is  that  it  is  largely  based  on  an  analyti¬ 
cal  derivation  of  the  Jacobian  matrix,  which  therefore  provides  exact  derivatives  of 
the  secondary  magnetic  held  with  respect  to  all  unknown  quantities  while  reducing 
the  computation  time. 

Inherent  to  the  method  itself,  various  initial  guesses  to  seed  the  algorithm  and 
maximum  iterations  numbers  have  been  tested,  depending  on  the  complexity  of  the 
cases  studied.  For  single  targets,  10  iterations  were  in  general  sufficient  whereas 
up  to  200  iterations  were  run  for  three-target  cases.  An  increased  number  of  it¬ 
erations  naturally  results  in  an  increased  computation  time,  which  reached  almost 
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0.1  1  10 

Time  [ms] 


(a)  Mortar  and  nose-cone. 


(b)  Baseplate. 

Figure  9:  Three-target  inversions  from  TEMTADS  data.  The  mortar  and  the  baseplate  are  held 
fixed  at  (jc,y,z)  =  (0,0,60)  [cm]  and  (jc,y,z)  —  (50,0,49)  [cm],  respectively,  while  the  nose-cone 
moves  along  x  from  case  to  case.  True  and  inverted  positions  are  given  in  Table  7. 
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3  minutes  in  case  of  TEMTADS  data  inverted  on  three  targets  (we  recall  that  the 
TEMTADS  data  set  is  typically  much  larger  than  the  MPV  data  set  so  that  inver¬ 
sions  are  more  time  consuming).  The  final  solution  produced  by  the  algorithm  is 
remarkably  insensitive  to  the  initial  guess  (apart  from  obvious  divergences),  pro¬ 
vided  that  the  step  at  each  iteration  is  weighted  by  a  parameter  0  <  a  <  1.  A  value 
of  a  =  0.5  yielded  very  consistent  inverted  positions  in  all  examples  shown  de¬ 
spite  large  variations  in  initial  guesses:  within  ±50  cm  for  both  x  and  y  coordinates 
and  between  —50  cm  and  —20  cm  for  the  depths.  The  inverted  polarizabilities, 
however,  were  more  sensitive  especially  in  the  three-target  cases,  yet  still  within 
acceptable  bounds.  To  resolve  possible  variations,  the  algorithm  was  run  on  typ¬ 
ically  10  realizations.  The  selection  of  the  final  solution  was  performed  by  using 
two  often  corroborating  criteria.  First,  the  histogram  of  solutions  for  the  inverted 
positions  was  examined  and  the  most  frequent  solution  was  declared  final.  Sec¬ 
ond,  the  error  was  computed  between  the  measured  held  and  the  computed  held 
and  the  solution  that  produced  the  lowest  error  was  declared  hnal.  In  all  cases  con¬ 
sidered  these  two  criteria  selected  identical  solutions.  It  should  also  be  mentioned 
that  simpler  cases  such  as  single  target  inversion  yielded  similar  hnal  solutions  in 
almost  all  the  cases,  so  that  the  criteria  proposed  above  were  mostly  useful  when 
inverting  multi-target  conhgurations. 

Finally,  we  emphasize  that  the  assumption  of  a  known  number  of  targets  (i.e. 
known  N  in  Eq.  (4.4))  is  necessary  for  the  development  of  the  algorithm  but  is 
not  required  in  practice  and  only  the  knowledge  of  an  upper  bound  is  needed.  If 
the  algorithm  is  run  with  N  larger  than  the  actual  number  of  targets,  the  algorithm 
may  yield  two  types  of  results.  First,  some  dipoles  cluster  around  a  certain  loca¬ 
tion  and  capture  the  secondary  held  of  a  given  UXO  better  than  a  single  dipole. 
This  situation  is  not  very  common  but  can  occur  for  very  heterogeneous  targets. 
More  commonly,  clusters  are  very  tight  and  the  dipoles  hnd  themselves  practically 
at  the  same  location.  The  corresponding  time-dependent  polarizabilities  need  to 
be  added  and  collectively  represent  the  signature  of  the  target  (individual  polariz¬ 
abilities  typically  do  not  exhibit  any  physical  behavior).  The  second  type  of  re¬ 
sults  obtained  if  N  is  larger  than  the  actual  number  of  targets  is  that  some  dipoles 
converge  to  apparently  irrelevant  locations,  with  polarizabilities  a  few  orders  of 
magnitude  lower  than  those  of  the  targets  themselves.  We  attribute  this  behavior 
to  the  algorithm’s  attempt  to  hnd  a  best  fit  for  noisy  data,  which  is  obtained  with 
a  few  well-identified  dipoles  and  other  weaker  ones  that  capture  average  noise 
contributions.  At  the  target  identification  stage,  these  spurious  dipoles  are  easily 
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filtered  out  due  to  their  unphysical  and  weak  responses.  Note  also  that  this  latter 
behavior  prevents  the  use  of  the  misfit  error  to  determine  the  number  of  targets 
since  smaller  errors  are  typically  obtained  with  more  dipoles.  Consequently,  the 
algorithm  does  not  require  the  exact  knowledge  of  the  number  of  targets  for  their 
proper  identification,  and  only  an  upper  bound  needs  to  be  estimated.  This  topic 
is  currently  under  investigation. 

4.3  Joint  Diagonalization  and  the  Orthonormalized  Volume  Magnetic  Source 


Methods 


We  apply  a  volume  dipole  distribution  model  to  next-generation  sensor  data  for 
multi-object  data  inversion  and  discrimination.  Discrimination  between  UXO  and 
harmless  objects  is  particularly  difficult  in  highly  contaminated  sites  where  two 
or  more  objects  are  simultaneously  present  in  the  held  of  view  of  the  sensor  and 
produce  overlapping  signals.  The  first  step  in  overcoming  this  problem  is  esti¬ 
mating  the  number  of  targets.  In  this  report  an  orthonormalized  volume  magnetic 
source  (ONVMS)  approach  is  introduced  for  estimating  the  number  of  targets, 
along  with  their  locations  and  orientations.  The  technique  is  based  on  the  discrete 
dipole  approximation,  which  distributes  dipoles  inside  the  computational  volume. 
First,  a  set  of  orthogonal  functions  are  constructed  using  fundamental  solutions 
of  the  Helmholtz  equation  (i.e.,  Green’s  functions).  Then,  the  scattered  magnetic 
held  is  approximated  as  a  summation  of  these  orthogonal  functions.  The  magni¬ 
tudes  of  the  expansion  coefficients  are  determined  directly  from  the  measurement 
data  without  solving  an  ill-posed  inverse-scattering  problem.  The  expansion  coef¬ 
ficients  are  then  used  to  determine  the  amplitudes  of  the  responding  volume  mag¬ 
netic  dipoles.  The  algorithm’s  superior  performance  and  applicability  to  live  UXO 
sites  are  illustrated  by  applying  it  to  the  bi-static  TEMTADS  multi-target  data  sets 
collected  by  NRL  personnel  at  the  Aberdeen  Proving  Ground  UXO  test- stand  site. 

4.3.a  Introduction 

Remediation  of  lands  contaminated  with  unexploded  ordnance  (UXO)  at  Depart¬ 
ment  of  Defense  and  Department  of  Energy  sites  has  been  identihed  as  one  of  the 
military’s  most  pressing  environmental  problems  [50].  As  a  result  of  past  military 
and  weapon-testing  activities,  UXO  are  found  at  both  active  and  formerly  used  de¬ 
fense  sites  (FUDS)  (closed,  transferred  and  transferring  ranges,  munitions  burning 
and  open  detonation  areas).  In  the  United  States  alone,  more  than  900  sites,  about 
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1 1  million  acres  of  land,  are  potentially  contaminated  with  UXO.  The  costs  of 
excavating  all  geophysical  anomalies  are  well-known  and  are  one  of  the  greatest 
impediments  to  the  efficient  cleanup  of  UXO,  particularly  at  highly  contaminated 
UXO  sites,  where  multiple  objects  are  often  present  simultaneously  in  the  sen¬ 
sors’  held  of  view.  State-of-the-art  electromagnetic  induction  (EMI)  sensors  have 
recently  emerged  as  the  most  cost-effective  sensing  technology  for  detecting  and 
discriminating  between  UXO  and  non-hazardous  items.  These  sensors  are  easy 
to  operate  in  various  environments  and  can  be  deployed  hand-held  or  on  a  variety 
of  mobile  platforms.  For  the  low  frequencies  and  high  conductivities  of  interest 
here,  induced  conduction  currents  are  much  stronger  than  displacement  currents, 
so  the  latter  can  usually  be  neglected  in  both  target  and  the  surrounding  soil.  In 
this  frequency  regime,  electric  fields  are  typically  quite  small,  and  the  conductiv¬ 
ity  of  the  surrounding  soil  is  typically  many  orders  of  magnitude  less  than  that  of 
the  metallic  target.  Together  with  the  low  frequency  range,  this  means  that  the 
dielectric  properties  of  the  surrounding  media  are  negligible  and  could  be  treated 
as  free  space,  which  makes  EMI  sensors  more  attractive  when  compared  to  elec¬ 
tromagnetic  wave  sensing  technologies. 

In  current  practice,  most  EMI  sensors  are  composed  of  separate  transmitting 
and  receiving  coils.  When  the  operator  activates  the  sensor,  a  current  runs  through 
the  transmitter  coils,  which  results  in  the  establishment  of  a  (“primary”  or  “princi¬ 
pal”)  magnetic  held  in  the  surrounding  space.  By  Faraday’s  law,  this  time-varying 
magnetic  held  induces  eddy  currents  in  highly  conducting  bodies;  ferromagnetic 
bodies  also  have  their  magnetization  affected  by  the  impinging  held.  These  cur¬ 
rents  and  magnetization  in  turn  generate  a  (“secondary”  or  “scattered”)  magnetic 
held  that  also  varies  with  time  and  induces  measurable  currents  in  the  receiving 
coils.  At  the  end,  the  electromagnetic  data  are  inverted  using  different  forward 
models. 

Recently,  several  EMI  sensing  and  data-inversion  techniques  [4,  20,  24,  50-64] 
have  been  developed  for  detecting  and  discriminating  between  UXO  and  non- 
UXO  items.  Typically  the  hrst  step  of  these  methods  is  the  recovery  of  a  set 
of  parameters  that  specify  a  physics-based  model  representing  the  object  under 
interrogation.  For  example,  in  EMI  sensing,  the  recovered  parameters  consist 
of  the  object’s  location  and  spatial  orientation  in  addition  to  “intrinsic”  param¬ 
eters  such  as  the  polarizability  tensor  (along  with  some  parametrization  of  its 
time-decay  curve)  in  dipole  models  or  the  amplitudes  of  responding  magnetic 
sources  in  the  NSMS  and  SEA  models.  EMI  responses  depend  nonlinearly  on 
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Tx 


Figure  10:  Metallic  objects  under  the  transmitter.  The  target’s  EMI  response  at  the  receiver  coil 
can  be  calculated  from  the  volume  magnetic  dipole  moment  dm. 


the  subsurface  object’s  location  and  orientation,  and  therefore  the  determination 
of  a  buried  object’s  orientation  and  location  is  a  non-linear  problem.  In  the  last 
decade  several  inverse  scattering  approaches  have  been  adapted  for  EMI  data  in¬ 
version  [52,  53,  55,  65,  66].  These  algorithms  are  carried  out  by  determining  an 
objective  function  as  a  goodness-of-ht  measure  between  the  forward  model  and 
measurement  data. 

Although  various  forward  EMI  models  (such  as  single-dipole,  multi-dipole  [20, 
50-53],  and  normalized  surface  magnetic  source  [4,  24])  have  demonstrated  great 
UXO  discrimination  capabilities  for  well-isolated  single  targets  using  data  pro¬ 
vided  by  next-generation  EMI  sensors,  single-target  discrimination  capability  is 
considerably  limited  for  overlapping  EMI  signals;  i.e.,  when  a  sensor  detects  EMI 
signals  from  two  or  more  targets  at  the  same  time.  The  obvious  extension  of 
current  EMI  modeling  to  accommodate  overlapping  EMI  signatures  is  multiple- 
target  modeling.  However,  when  we  apply  multi-target  modeling  to  real  data,  the 
accuracy  of  the  solution  is  generally  limited  by  the  ill-posed  nature  of  the  inverse 
problem  and,  thus,  multi-target  modeling  with  more  than  two  targets  is  seldom 
applied  in  practice  in  UXO  detection  and  discrimination. 

Even  when  multi-target  EMI  forward  modeling  is  applied  to  high-fidelity  next- 
generation  sensor  EMI  data,  the  first  step  is  estimating  the  number  of  targets  and 
then  extracting  the  intrinsic  parameters  of  these  multiple  targets,  which  generally 
requires  a  highly  multidimensional  linear  and  nonlinear  optimization  search.  In 
many  cases,  extracting  linear  parameters  from  a  limited  set  of  data  is  an  ill-posed 
problem  with  many  solutions  and  for  which  no  existing  technique  can  guarantee 
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that  the  true  solution  will  be  attained  within  the  practical  limits  of  computational 
time.  Therefore,  the  success  of  multiple-target  parameter  estimation  greatly  de¬ 
pends  on  how  closely  a  forward  model  with  a  small  number  of  model  parame¬ 
ters  can  take  into  account  the  underlying  physics  of  EMI  scattering  phenomena. 
One  way  to  incorporate  the  physics  faithfully  is  to  apply  the  orthonormalized  vol¬ 
ume  magnetic  source  (ONVMS)  model  to  the  multi-object  UXO  discrimination 
problem.  The  algorithm  is  based  on  the  discrete  interacting  dipole  approxima¬ 
tion  and  represents  EMI  scattered  fields  using  a  set  of  orthogonal  functions.  This 
paper  combines  the  ONVMS  with  a  global  search  technique  called  differential 
evolution  [65,  66]  to  determine  target  locations  and  orientation  starting  from  next- 
generation  sensor  data.  It  then  uses  joint  diagonalization  [67]  to  estimate  the  num¬ 
ber  of  contributing  targets  and  to  de-noise  signals. 

The  paper  is  organized  as  follows:  Sec.  4.3.b  presents  the  orthonormalized  vol¬ 
ume  magnetic  source  technique,  Sec.  4.3.C  describes  orthonormalized  functions, 
and  Sec.  4.3.d  presents  the  joint  diagonalization  approach.  Section  4.3.e  shows 
experimental  and  numerical  results  that  demonstrate  the  applicability  of  the  com¬ 
bined  ONVMS -DE  algorithm  and  joint  diagonalization  approach  to  the  detection 
and  discrimination  of  multiple  buried  objects  starting  from  actual  EMI  data  col¬ 
lected  by  the  TEMTADS  array. 

4.3.b  The  Orthonormalized  Volume  Magnetic  Source  Mode 

Let  us  consider  permeable  and  highly  conducting  metallic  3-D  objects  placed  in 
a  time  varying  primary  magnetic  held  (see  Fig.  10).  As  is  well  established  in  the 
EMI  frequency  regime,  we  assume  that  displacement  currents  dD/dt  are  negligi¬ 
ble.  Further,  the  conductivity  of  the  soil  and  the  electric  fields  within  it  are  small 
and  combine  to  produce  negligible  currents.  The  primary  magnetic  held  pene¬ 
trates  the  object  to  some  degree,  inducing  eddy  currents/magnetic-dipoles  inside 
it.  In  return,  these  induced  sources  produce  a  secondary  or  scattered  magnetic 
held  outside  the  objects  that  can  be  represented  as  a  sum  of  magnetic  helds  due  to 
a  normalized  set  of  magnetic  sources  distributed  over  a  volume: 


(4.14) 


v 


A 

where  p  =  t,f  is  time  or  frequency,  R  is  the  unit  vector  along  R  =  r  —  r  v,  r  v  is 
the  position  of  the  v'-th  inhnitesimal  dipole  in  the  volume  V,  r  is  the  observation 
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point  (see  Fig.  10),  and  I  is  the  identity  dyad.  In  the  ONVMS  approach  the  in¬ 
duced  magnetic  dipole  moment  m(rv/,  p)  at  point  ly  on  the  surface  is  related  to  the 
primary  held  via  m(rv/,p)  =  M(rv/,p)  •  Hpr(r v/),  where  M(rv/,p)  is  the  symmet¬ 
ric  polarizability  tensor.  Overall,  in  the  NVMS  approach  the  scattered  magnetic 
held  is  approximated  as  a  superposition  of  the  helds  radiated  by  these  elementary 
sources,  whose  amplitudes  are  normalized  with  the  impinging  primary  magnetic 
held  H/;r.  According  to  a  simplihcation  of  the  elementary  atomic  model  of  matter, 
all  materials  are  composed  of  atoms,  each  with  a  positively  charged  nucleus  and  a 
number  of  orbiting  negatively  charged  electrons.  The  orbiting  electrons  cause  cir¬ 
culating  currents  and  form  microscopic  magnetic  dipoles  inside  the  targets.  In  the 
absence  of  an  external  magnetic  held  the  magnetic  dipoles  of  atoms  of  most  mate¬ 
rials  have  random  orientations,  resulting  in  no  magnetic  moment.  The  application 
of  an  external  time-varying  magnetic  held,  by  Faraday’s  law,  induces  eddy  cur¬ 
rents  in  highly  conducting  bodies  by  an  alignment  of  the  magnetic  moments  of  the 
spinning  electrons  and  a  magnetic  moment  due  to  a  change  in  the  orbital  motion 
of  electrons.  These  currents  and  magnetization  in  turn  generate  a  (“secondary” 
or  “scattered”)  magnetic  held  that  also  varies  with  time  and  induces  measurable 
currents  in  the  receiving  coils.  Since  the  induced  magnetic  dipoles/eddy  currents 
are  distributed  inside  the  object,  the  spatial  distribution  of  the  magnitudes  of  the 
dipoles  indicate  the  locations  and  orientations  of  any  targets  present  inside  the 
volume.  The  equation  (4.14)  in  discrete  form  can  be  written  as 

Nv  i  _  Nv 

H(r)  -  V  - - -J (3R,R;  - 1)  •  mt  =  Y  G,(r)  ■  mu  (4. 1 5) 

i=  i  47rMo ,■= t 
where  G;( r)  is  Green’s  dyad. 


4.3.C  Orthonormalized  Functions 

The  scattered  magnetic  held  at  any  point  outside  the  volume  (4.14)  can  be  repre¬ 
sented  using  orthonormalized  functions  as 

Nv 

H(r)  =  £i URi)K  (4.16) 

i=  1 

where  R/  =  r  -  r J,  {rj }^=1  £  V,  r)  is  an  orthonormalized  function,  and  bj  is 
orthogonal-function  expansion  coefficient.  The  orthonormalized  functions  ij//(R;) 
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are  linearly  independent  at  the  observation  point: 

M 

£fec„  =  0.  (4.17) 

n=  1 

This  means  that  all  coefficients  c\  =  C2  =  . . .  =  cm  =  0.  In  addition,  the  set  of 
orthogonal  functions  is  complete,  which  guarantees  that  there  is  a  set  of  coeffi¬ 
cients  cn,  n  =  1, . . .  ,M,  that  will  specify  the  magnetic  held  H  at  any  point  with 
arbitrary-eg  accuracy: 


N 


M 


n—  1 


2 

ds  <  £1where£  <C  1. 


The  orthonormal  functions  1/4  are  constructed  as  linear  combinations  of  the  Green’s 
functions  using  the  Gram-Schmidt  orthonormalization  process, 


m—  1 

Vm  Gm  % '  A-mki 

k=l 


(4.18) 


where  Gm  is  Green’s  dyad  from  (4.15)  and  the  tensor  A, „^is  determined  as 


n—\ 


A-mn  —  G  |  Cnm  A-nk  '  F k  '  ^mk  )  j  i^-mn  —  O5  m  — 

fc=l 


(4.19) 


The  functions  i/4satisfy  the  orthogonality  conditions 

[  (fl-fk)  dv  = 


0,  m  7^  k 
Fm,  m  =  k 


and 


Gm  '  Gk  )  dv  Cmk- 


(4.20) 


(4.21) 


The  orthogonal  functions  expansion  coefficients  bj  in  equation  (4.16)  can  be  de¬ 
termined  from  actual  data  Hdata  using  the  orthogonality  property  of  the  ij/  given 
in  (4.20): 

bi  =  Pr1  ■  I  iff  •  Hdata  dv.  (4.22) 


In  addition,  (4.18)  can  be  rewritten  as 


Vi(s)  —  ^  Gn (s)  • Bin i 
11—  1 


(4.23) 
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where  the  tensor  coefficient  Bm  can  be  determined  using  the  tensors  Am-  from 
(4.19)  through 

B  m-  Y  Anj  •  B,nn ,  Bmm  =  /,  =  Am(m_iy  (4.24) 

n=  1 

m  >  2 
i  <  m 

Finally,  by  combining  (4.16)  and  (4.23),  the  scattered  magnetic  held  can  be  ex¬ 
pressed  as 

Nv  Nv  i  _  i  /  m  \  Nv 

H(r)  =  £  pi(s)  •*,•  =  ££  <5„(r)  •  4  •  7>;  =  £  <3,(r)  •  £  4  •  b„  =  £  G,(r)  • 

/=1  i=ln=l  /=!  \/i=l  /  /=1 

(4.25) 

and  comparing  (4.15)  to  (4.25)  and  using  the  linear  independence  of  Green’s  func¬ 
tions  one  can  conclude  that 

i 

Bin  ‘  bi . 

72  =  1 

Thus  the  amplitudes  of  the  responding  magnetic  dipoles  are  determined  without 
solving  a  linear  system  of  equations. 


4.3.d  Joint  Diagonalization 


Recently,  under  the  Strategic  Environmental  Research  and  Development  Program 
(SERDP),  next  generation  EMI  sensors  have  been  developed  for  discriminating 
between  subsurface  UXO  and  clutter.  One  such  sensor  is  the  Time  domain  Electro¬ 
magnetic  Multi-Sensor  Towed  Array  Detection  Systems  (TEMTADS).The  system 
consists  of  25  transmit/receive  pairs  of  square  coil  antennas  arranged  in  a  5  x  5 
grid,  each  consisting  of  a  35-cm  transmitter  loop  and  a  25-cm  receiver  loop.  The 
sensor  activates  the  transmitter  loops  in  sequence,  one  at  a  time,  and  for  each  trans¬ 
mitter  all  receivers  receive,  measuring  the  complete  transient  response  over  a  wide 
dynamic  range  of  time:  from  approximately  100  /is  to  25  ms  distributed  in  123 
time  gates.  The  sensor  thus  provides  625  spatial  data  points  at  each  instrument 
location,  offering  unprecedented  positional  accuracy.  Recent  TEMTADS  blind 
discrimination  studies  [24]  show  perfect  performance  when  the  number  of  targets 
was  known  a  priori.  However,  determining  the  number  of  potential  targets  and 
locating  them  is  still  a  challenging  problem.  To  take  advantage  of  the  high  fidelity 
data  provided  by  TEMTADS,  here  we  employ  the  joint  diagonalization  algorithm 
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to  estimate  the  number  of  targets  directly  from  the  data  using  eigenvalues.  Let  us 
briefly  describe  the  joint  diagonalization  algorithm  for  the  TEMTADS  system.  For 
a  given  time  channel  tq,  q  =  1,2,...,  Nq,  where  Nq  is  the  number  of  time  channels, 
let  us  define  H^m  as  the  z-component  of  the  magnetic  held  measured  by  the  m-th 
receiver  coil  when  the  k-th  transmitter  is  active,  and  the  matrix 


fi(k)  = 


#  1,2  • 

••  #i,m 

#2,1 

#2,2  • 

■  ■  #2,M 

Hr,  i 

Hk,  2  ' 

•  •  Hk.m 

(4.26) 


as  a  set  of  measured  data/vector  for  the  k-th  transmitter,  where  k  =  1,  2,  . . . ,  K. 


Figure  11:  A  schematic  diagram  of  the  three-target  TEMTADS  experimental  setup.  The  targets 
are:  1)  4.2”  mortar  oriented  horizontally  along  the  x  axis  and  placed  at  (x=y=0)  60  cm  bellow  the 
TEMTADS  centerline;  2)  A  half-shell  oriented  horizontally,  placed  at  x=-50  cm,  y=0,  and  47  cm 
bellow  the  TEMTADS;  3)  A  base  plate  oriented  vertically,  placed  at  x=y=0,  and  27  cm  bellow  the 
TEMTADS.  The  TEMTADS  center  is  at  x=y=0;  z=60  cm. 


For  the  TEMTADS  system  M  -  25  is  the  number  of  receivers  and  K  =■  25  is  the 
number  of  transmitters.  For  each  time  channel  the  M  x  K  matrix  [HcI\  can  be 
expressed  with  eigenvectors  U (tq)  and  eigenvalues  D{tq)  using  the  singular  value 
decomposition: 

fi(tq)  =  0(tq)b(tq)VT(tq).  (4.27) 

However,  in  order  to  determine  the  number  of  potential  targets,  we  need  to  find  an 
eigenvector  V  that  will  be  shared  by  all  ^Hd{tq)  j  ;  ^matrices  and  will  also  make 
vanish  all  off-diagonal  elements  of  the  matrix 
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D(t,)  =  PTti(tq)V,  q  =  1  to  Nq. 


(4.28) 


In  general,  it  is  not  the  case  that  V  will  cancel  all  off-diagonal  elements  of 


D(tq).  The  technique  that  finds  a  unitary  V  matrix  that  minimizes  the  Hd(t, 


w, 


q=  i 


matrices’  off-diagonal  elements  is  called  the  joint  diagonalization  approach  [67]. 


10'° 

Time  [sec] 


10 


b) 


Figure  12:  Eigenvalues  versus  time  for  4.2”  mortar  TEMTADS  data:  a)  All  25  eigenvalues,  b)  The 
three  highest  eigenvalues  are  associated  with  the  4.2”  mortar  and  the  next  three  correspond  to  non 
dipole  terms. 


Time  [sec]  Time  [sec] 


Figure  13:  Eigenvalues  versus  time  for  three-target  TEMTADS  data:  b)  All  25  eigenvalues,  b)  The 
first  three  highest  eigenvalues  are  associated  with  4.2”  mortar,  half  shell  and  base  plates. 
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4.3.e  Results 

First  we  apply  the  joint  diagonalization  model  to  measured  4.2”  mortar  TEM- 
TADS  data  collected  over  a  calibration  test  stand.  The  object  was  placed  hori¬ 
zontally  at  a  60  cm  depth  under  the  TEMTADS  array  centerline.  The  H(tq)  data 
were  diagonalized  using  (4.28).  Figure  12(a)  shows  all  25  eigenvalues  vs.  time 
for  the  TEMTADS  magnetic  held  H{tq).  The  result  shows  that  there  are  three 
strong  eigenvalues  and  that  the  rest  are  small.  Since  the  data  was  taken  only  for 
one  4.2”  mortar,  we  could  conclude  that  these  strong  eigenvalues  correspond  to 
the  mortar  and  rest  of  the  eigenvalues  are  associated  to  sensor  noise.  Figure  12(b) 
shows  the  eigenvalues  larger  than  a  particular  threshold  after  smaller  eigenvalues 
are  removed. 


w 


o 

1 

H 


Time  [m  sec] 


Figure  14:  The  total  ONVMS  versus  time  for:  a)  4.2”  mortar  and  b)  half-shell  inverted  from 
three-object  TEMTADS  data,  when  base  plates  were  placed  at  the  six  different  locations  shown  in 
Table  15  below. 


Next  we  applied  the  proposed  approach  to  three-object  TEMTADS  data.  The 
schematic  diagram  of  the  experimental  setup  is  depicted  on  Fig.  11.  The  data  were 
acquired  by  NRL  personnel  at  the  Blossom  Point  UXO  test  stand  site.  The  sensor 
was  stationary,  and  the  4.2”  mortar  was  placed  horizontally  at  61  cm  under  the 
array  centerline  while  the  half  shell  was  lying  horizontally  and  placed  under  the 
array  centerline  at  47  cm  depth,  x  =  -50  cm  offset,  y  =  0.  The  base  plate  was  kept  at 
27  cm  depth  and  moved  from  the  center  along  negative  x-axis  at  10  cm  increments 
while  keeping  y  =  0.  First  we  used  the  joint  diagonalization  technique  to  estimate 
the  number  of  potential  targets.  Figure  13(a)  shows  all  25  eigenvalues  vs.  time  for 


10 

Time  [m  sec] 
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Case  # 

Base  plate  location 
inverted/  actual 

[cm] 

X 

z 

1 

-5/0 

-27/-30 

2 

-12/-10 

-34/-30 

3 

-22/-20 

-35/-30 

4 

-35/-30 

-37/-30 

5 

-42/-40 

-37/-30 

6 

-54/-50 

-36/-30 

Figure  15:  Inverted  and  actual  locations  of  the  base  plate. 


Time  [m  sec] 


Figure  16:  Inverted  total  ONVMS  versus  time  for  the  baseplate  at  six  different  locations  along  the 
x-axis. 
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Time  [sec]  Time  [sec] 


Figure  17:  Eigenvalues  versus  time  for  TEMTADS  data  corresponding  to  three  targets,  a  105-mm 
projectile  and  two  spheroids:  a)  All  25  eigenvalues,  b)  Nine  highest  eigenvalues. 


the  TEMTADS  magnetic  field  H(tq)  when  the  base  plate  was  placed  at  x  =  y  =  0 
and  the  depth  was  27  cm.  The  result  shows  that  there  are  nine  strong  eigenvalues 
corresponding  to  the  mortar,  half  shell,  and  base  plates  (see  Figure  13(b)).  The 
rest  of  the  eigenvalues  are  small  and  are  again  associated  to  sensor  noise,  just  as 
in  Figure  12.  We  removed  the  sensor  noise  from  the  data  by  removing  the  small 
eigenvalues  and  their  corresponding  eigenvectors.  Then  we  applied  the  combined 
ONVMS  and  differential  evolution  algorithm  [65,  66]  to  the  de-noised  data  to 
determine  target  locations  and  the  corresponding  total  ONVMS  in  order  to  perform 
discrimination.  The  inverted  locations  of  the  base  plates  are  summarized  in  Table 
1.  The  results  show  that  the  algorithm  inverted  the  small  object’s  burial  depth 
and  location  with  5-cm  accuracy  for  all  six  cases.  The  combined  ONVMS-DE 
algorithm  also  inverts  the  locations  of  the  4.2”  mortar  and  half-shell  with  3  cm 
accuracy  for  all  six  case  (the  results  are  not  shown  here). 

The  algorithm  outputs  the  total  ONVMS  versus  time  for  each  targets.  Fig¬ 
ure  14  shows  the  total  ONVMS  for  4.2”  and  half-shell  for  all  six  cases.  The  results 
demonstrate  that  the  algorithm  extracts  a  stable  total  ONVMS  decay  curves  for  all 
three  targets  for  the  six  different  configurations;  see  Figure  14  and  Figure  16.  Thus 
the  algorithm  is  capable  to  classify  correctly  between  UXO  and  non-UXO  items. 

Finally,  to  demonstrate  the  applicability  of  the  joint  diagonalization  and  ONVMS- 
DE  algorithm  for  a  highly  contaminated  UXO  scenario  in  which  a  big  target  dom¬ 
inates  over  small  targets,  again  we  employ  the  method  on  a  three-target  TEM¬ 
TADS  data  set.  The  targets  are  a  105-mm  projectile  of  length  49  cm  and  two 
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Figure  18:  Inverted  total  ONVMS  versus  time  for  the  105  mm  projectile. 


steel  spheroids  with  major  to  minor  axis  ratio  bla  =  4  and  a  =  2.5cm.  The  sensor 
was  stationary,  and  the  105  mm  projectile  was  placed  horizontally  at  59  cm  under 
the  array  centerline,  while  the  spheroids  were  placed  at  depth  33  cm  and  x  =  0, 
y  =  10  cm  and  x  =  0y  =  -10cm  respectively.  We  apply  the  joint  diagonalization 
technique  to  estimate  the  number  of  potential  targets.  Figure  17(a)  shows  all  25 
eigenvalues  vs.  time  for  the  TEMTADS  magnetic  field.  The  result  shows  that 
even  when  the  spheroids  are  close  to  the  sensor  their  corresponding  eigenvalues 
are  orders  of  magnitudes  less  than  those  of  the  105-mm  projectile.  This  indicates 
that  the  contribution  from  the  spheroids  is  negligible  (Fig.  17(b)).  We  removed 
the  sensor  noise  and  the  spheroids’  signals  from  the  data  by  making  all  small 
eigenvalues  and  their  corresponding  eigenvectors  be  zero.  Again  we  applied  the 
combined  ONVMS/differential  evolution  algorithm  to  the  remaining,  de-noised 
data  for  determining  the  target’s  location  and  the  corresponding  total  ONVMS  for 
discrimination.  The  algorithm  correctly  predicted  the  location.  The  approach  in¬ 
verted  a  stable,  symmetric  total  ONVMS  versus  time  for  the  105-mm  projectile, 
as  shows  Fig.  18. 
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4.3.f  Conclusion 

The  orthonormalized  volume  magnetic  source  technique  was  introduced  for  multi¬ 
object  detection  and  discrimination.  The  algorithm  was  applied  to  next-generation 
EMI  array  sensors,  in  particular  to  the  TEMTADS  array.  The  ONVMS  was  com¬ 
bined  with  DE  to  estimate  target  depths  and  locations.  In  addition,  a  joint  diago- 
nalization  algorithm  was  adapted  to  TEMTADS  data  for  estimating  the  number  of 
potential  targets.  The  TEMTADS  sensor,  a  bi-static  system  with  25  x  25  Tx/Rx 
combinations,  yields  a  25  x  25  matrix  for  each  time  channel.  First,  by  examining 
the  eigenvalue  structure  of  this  matrix  we  can  efficiently  estimate  the  number  of 
targets,  and  then  by  utilizing  the  combined  ONVMS-DE  algorithm  we  can  locate 
the  subsurface  targets  and  efficiently  extract  their  intrinsic  parameters  (such  as  the 
polarizability)  for  discrimination  purposes.  The  accuracy  of  the  combined  tech¬ 
nique  was  tested  against  time-domain  TEMTADS  data  collected  at  calibration  and 
test  stand  sites.  The  results  show  that  that  the  combined  model  accurately  predicts 
the  number  of  targets,  their  burial  depths  and  locations,  andstable  discrimination 
feature  parameters  that  could  be  used  to  discriminate  between  objects  of  interest 
and  non-hazardous  targets. 

5  Data,  Results,  and  Discussion 

5.1  TEMTADS  Data  and  Results 

The  spheroids  detailed  in  Sec.  4.1  are  now  all  combined  and  in  the  custody  of 
Dan  Steinhurst  at  NOVA  Research  who  maintains  and  operates  the  TEMTADS 
instrument  under  a  contract  from  NRL.  As  part  of  this  project,  an  extensive  set  of 
data  was  acquired  over  many  configurations  of  multiple  targets,  including  these 
canonical  targets.  Table  5.1  gives  a  sample  of  the  data  collected.  For  a  full  list, 
contact  Ben  Barrowes. 

Data  cases  1 17-126  were  given  to  Dr.  Fridon  Shubitidze  and  Dr.  Tomasz  Grze- 
gorczyk  as  a  blind  data  set.  Figure  19  shows  a  sample  picture  of  the  TEMTADS 
and  the  targets.  Using  the  JD  and  ONVMS  methods  (see  Sec.  4.3),  he  was  able  to 
correctly  identify  all  targets  in  this  blind  data  set  (see  Fig.  20). 
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5.2  MetalMapper  Data  and  Results 


Figure  19:  Sample  setup  for  TEMTADS  multitarget  data  collection. 


5.2  MetalMapper  Data  and  Results 


To  date,  there  is  no  good  set  of  data  from  the  MetalMapper  over  multiple  targets. 
As  a  substitute,  we  constructed  a  set  of  multitarget  data  from  single  target  data 
via  superposition.  This  data  comes  from  test  data  that  Dave  George  took  while 
calibrating  the  MetalMapper.  For  this  test,  one  of  three  targets: 

1.  10cm  diam  sphere 


2.  13cm  piece  of  rebar 

3.  21cm  diameter  copper  ring 


were  placed  in  one  of  five  positions  under  the  MetalMapper  (see  Fig.  21). 

In  the  magnetoquasistatic  (MQS)  regime,  the  data  can  be  superposed  without 
difficulty.  Therefore,  multitarget  test  data  was  synthesized  by  adding  up  the  data 
from  different  cases  of  targets  in  different  quadrants.  Table  9  shows  the  original 
data  from  which  these  blind  tests  were  synthesized.  These  original  data  were 
superposed  (added)  in  the  following  combinations  (see  Table  10)  before  being 
given  to  Dr.  Fridon  Shubitidze  and  Dr.  Tomasz  Grzegorczyk  as  a  blind  data  set. 

Using  the  multiple  dipole  Jacobian  based  gradient  search  method  (see  Sec.  4.2), 
Dr.  Grzegorczyk  was  able  to  correctly  identify  all  the  targets  in  the  semiblind  data 
sets  (see  Figs.  22  and  23).  The  results  for  each  case  are  computed  in  under  one 
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5.3  MPV2  Data 


Case  # 

Trg#l 

Trg#2 

Trg-3 

Trg#4 

Trg*5 

Trg#6 

117 

60mm 
(0, -40,-3  8) 

X 

X 

X 

X 

X 

118 

60mm 
(0,-40, -38) 

SSI 

(21,22,-26) 

X 

X 

X 

X 

119 

60mm 
(0, -38.-38) 

SSI 

(20,20,-26) 

SP4S 

(-20,-20,-28) 

X 

X 

X 

120 

60mm 
(0, -37,-40) 

SSI 

(22,21,-25) 

SP4S 

(-19,-19,-29) 

SP4A 

(-14,28,-30) 

X 

X 

121 

60mm 
(0,-34, -32) 

SSI 

(18,23,-26) 

SP4S 

(-19,-20,-31) 

SP4A 

(-24,19,-33) 

SP4S 

(30.-8.-28) 

122 

60mm 
(0.-35, -37) 

SSI 

(23,21,-21) 

SP4S 

(-10,-17,-33) 

SP4A 

(-24,-19,-26) 

SP4S 

(30,-8. -28) 

105  mm 
(0,37,  -57) 

123 

60mm 
(0,-35, -37) 

SSI 

(23,24,-21) 

SP4S 

(-12,-10,-33) 

SP4A 

(-28.-30,-24) 

X 

105  mm 
(0,41,  -57) 

124 

60mm 
(0.-37, -38) 

SSI 

(20.19,-19) 

SP4S 

(-19,-20,-29) 

X 

X 

105  mm 
(0,36,  -57) 

125 

60mm 
(0, -38,-39) 

SSI 

(22,20,-23) 

X 

X 

105  mm 
(0,38,  -58) 

126 

60mm 
(0, -38,-39) 

X 

X 

X 

105  mm 
(0,41,  -60) 

Figure  20:  Blind  inversion  results  from  Dr.  Shubitidze  using  the  combined  JD  and  ONVMS  meth¬ 
ods.  All  targets  correctly  identified  and  positions  are  within  the  size  of  the  target. 


minute  if  the  number  of  targets,  N,  is  known.  Finding  N  for  this  method,  however, 
still  requires  an  expert’s  interpretation  of  the  polarizabilities  assuming  different  N. 

Using  the  combined  JD  and  ONVMS  methods,  each  of  the  1-5  targets  was  also 
correctly  identified  for  this  blind  set  of  MetalMapper  data  (see  Fig.  24).  These 
results  take  longer  (10’s  of  minutes  for  each  case)  to  compute,  but  are  aided  by  an 
estimate  for  N  provided  by  the  JD  method. 

5.3  MPV2  Data 

MPV2  data  was  collected  recently  at  CRREL  in  Hanover,  NH  in  support  of  several 
SERDP  project  and  with  many  goals  in  mind: 

•  Collect  data  over  small  targets  at  several  depths  including  greater  than  11 
times  the  diameter  (in  support  of  MR-2106) 

•  Collect  data  over  2  and  3  target  scenarios  both  small  and  larger  targets  in 
several  combinations  (in  support  of  MR-2106  and  MR- 1664) 

•  Collect  data  over  deep,  larger  targets 
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5.3  MPV2  Data 


Figure  21:  MetalMapper  test  setup  with  one  of  three  targets  in  one  of  five  positions.  Blind  multi¬ 
target  test  data  synthesized  by  superposing  different  data  cases. 


•  Collect  data  of  small  and  large  targets  with  a  clutter  layer  between  the  targets 
and  the  MPV2  (in  support  of  MR-2106  and  MR-1664  and  MR-1005) 

•  To  carefully  characterize  the  beacon  positioning  system  accuracy  in  terms  of 
both  position  and  orientation 

•  To  investigate  the  integration  time  of  the  DAQ  window  of  the  MPV2  in  order 
to  optimize  the  data  acquisition  protocol  (in  support  of  MR- 1005  and  future 
MPV2  projects) 

•  To  collect  some  data  used  for  blind  tests  for  our  recent  algorithms  as  applied 
to  handheld  sensors 
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File# 

target 

depth  to 
bottom  of  target 

orientation 

quadrant 

2 

sphere 

35cm 

N/A 

1 

3 

ring 

35cm 

horizontal(H) 

1 

4 

ring 

35cm 

vertical(V)-  Odeg 

1 

5 

rebar 

35cm 

horizontal  -  Odeg 

1 

6 

sphere 

35cm 

N/A 

2 

7 

ring 

35cm 

H 

2 

8 

ring 

35cm 

V  -  Odeg 

2 

9 

rebar 

35cm 

H  -  Odeg 

2 

10 

sphere 

35cm 

N/A 

3 

11 

ring 

35cm 

V  -  Odeg 

3 

12 

ring 

35cm 

H  -  Odeg 

3 

13 

rebar 

35cm 

FI  -  Odeg 

3 

14 

sphere 

35cm 

N/A 

4 

15 

ring 

35cm 

H  -  Odeg 

4 

16 

ring 

35cm 

V  -  Odeg 

4 

17 

rebar 

35cm 

H  -  Odeg 

4 

18 

sphere 

35cm 

N/A 

0 

19 

ring 

35cm 

H  -  Odeg 

0 

20 

ring 

35cm 

V  -  Odeg 

0 

21 

rebar 

35cm 

H  -  Odeg 

0 

Table  9:  Original  MetalMapper  calibration  data. 


We  collected  this  data  in  February  2011,  and  we  believe  our  data  set  will  fulfill 
all  these  stated  goals  above,  but  we  have  not  had  time  to  produce  any  results  from 
this  data  set. 

6  Conclusions  to  Date 

SERDP  project  MR- 1664  entitled  “Isolating  and  Discriminating  Overlapping  Sig¬ 
natures  in  Cluttered  Environments”  is  approximately  halfway  complete.  We  have 
developed  two  major  methods  of  Ending  multiple  targets  in  EMI  data  even  when 
the  targets  are  in  close  proximity.  The  multiple  dipole  Jacobian  based  gradient 
search  method  (see  Sec.  4.2)  finds  multiple  targets  by  assuming  that  each  target 
exhibits  a  dipolar  response.  A  gradient  search  algorithm  then  searches  for  these 
multiple  sources.  At  this  point,  this  method  needs  to  know  the  number  of  targets  in 
order  to  perform  accurately.  The  JD  method  provides  this  estimate  for  N  to  either 
the  Jacobian  search  method  or  the  ON  VMS  method.  The  ONVMS  method  (see 
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blind  test  # 

original  file  #’s 

1 

[18,2] 

2 

[13,2,7] 

3 

[4,7,13,14] 

4 

[9,21] 

5 

[16,3] 

6 

[19,10,6,5] 

7 

[8,21,17] 

8 

[15,11,9] 

9 

[5,8,10,15,20] 

Table  10:  MetalMapper  blind  test  data  configurations  synthesized  by  superposing  data  cases  from 
Table  9. 


Case  nb 

Hb  of  targets 

Targets  id 

1 

2 

2 

18 

- 

- 

2 

3 

2 

7 

13 

- 

3 

3 

4 

7 

13 

14 

4 

2 

21 

9 

- 

- 

5 

2 

3 

16 

* 

- 

6 

4 

5 

6 

10 

19 

7 

3 

8 

21 

17 

- 

8 

3 

11 

9 

15 

- 

9 

5 

5 

8 

10 

15  20 

Figure  22:  Results  from  using  the  multiple  dipole  Jacobian  based  gradient  search  method.  Refer 
to  Tables  9  and  10. 


Sec.  4.3)  searches  for  dipole  like  (but  orthonormal)  sources  in  the  data  without 
inversion  difficulties  arising  from  solving  linear  systems.  Both  of  these  methods 
are  able  to  distinguish  at  least  six  different  targets  in  EMI  data  from  either  the 
TEMTADS  or  MetalMapper  instruments.  The  combined  JD  and  ONVMS  method 
(together  with  expert  guidance)  performed  very  well  at  the  Camp  Butner  demon¬ 
stration  site  at  the  Partners  2010  symposium.  We  still  plan  to  investigate  Particle 
Filtering  techniques  as  part  of  the  remainder  of  this  project. 
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MS*  MS* 

MalalmappM  llm—  (•  | 


Figure  23:  Polarizabilities  extracted  for  bind  case  9  using  the  Jacobian  multiple  dipole  search 
method.  Refer  to  Tables  9  and  10. 
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Case  # 

Object  1/ 
Location  (cm) 

Object  2/ 

Location  (cm) 

Object  3/ 

Location  (cm) 

Object  4/ 
Location  (cm) 

Object  5/ 
Location  (cm) 

1 

Sphere  / 

Sphere  / 

N/A 

N/A 

N/A 

(26, 28,  -37) 

(0,0,  -35) 

2 

Sphere  / 

Ring  / 

Rebar/ 

N/A 

N/A 

(26, 28,  -37) 

(-26, 21,  -28) 

(-27,  -28,  -38) 

3 

Ring  / 

Ring/ 

Rebar  / 

Sphere/ 

N/A 

(31, 28,  -29) 

(-26, 21,  -42) 

(-27,  -28,-38) 

(29,  -26,-35) 

4 

Rebar/ 

Rebar/ 

N/A 

N/A 

N/A 

(-28, 24,  -38) 

(0,0, -36) 

5 

Ring/ 

Ring/ 

N/A 

N/A 

N/A 

(26, 28,  -41.5) 

(29,-27,-29) 

6 

Rebar/ 

Sphere/ 

Sphere/ 

Ring/ 

N/A 

(28,23,-38) 

(-28, 24,-35) 

(-26,  -30,  -35) 

(0,0, -40) 

7 

Ring/ 

Rebar/ 

Rebar/ 

N/A 

N/A 

(-30,23,-28) 

(25,-32,37) 

(0,0, -36) 

8 

Rebar/ 

Ring/ 

Ring/ 

N/A 

N/A 

(-28,24,-38) 

(-25,-30,-28.5) 

(29.-27,-43) 

9 

Rebar/ 

Ring/ 

Sphere/ 

Ring/ 

Ring 

(28,  23,  -38) 

(-30,23,  -28) 

(-26,-30,  -30) 

(29,-27,-29) 

(0,0, -27) 

Figure  24:  Inversion  results  from  using  the  combined  JD  and  ONVMS  methods.  Refer  to  Tables  9 
and  10. 


Benjamin  Barrowes,  Fridon  Shubitidze 


-53- 


engineering 

AT  DARTMOUTH 


MR- 1664  -  March  2011  Interim  Report 


A  APPENDICES  -  SUPPORTING  DATA 


A  Appendices  -  Supporting  Data 

A.l  MPV  Single  Target  Blind  Test  Results 

We  collected  data  using  the  MPV  1  instrument  over  several  single  target  cases  as 
a  blind  test  for  our  models.  During  the  data  collection,  we  employed  the  laser 
positioning  system  to  locate  our  data  though  we  used  a  grid  as  a  rough  guide  when 
we  moved  the  MPV  instrument  itself.  We  took  data  at  two  heights,  separated  by  10 
cm  or  two  blue  board  thicknesses.  At  each  height  and  after  acquiring  background 
data,  we  took  a  total  of  17  measurement  in  a  four  by  four  grid  followed  by  a  single 
shot  located  precisely  on  our  grid.  We  used  this  last  data  point  to  define  a  global 
coordinate  system  for  each  set  of  17  data  points.  The  true  locations  and  orientation 
of  the  five  targets  are  given  in  Table  1 1 .  These  results  are  very  encouraging  for 
single  target  cases  even  considering  only  5  data  shots  (=75  data  points)  from  the 
MPV.  Sample  results,  presented  at  Partners  2007,  are  given  in  Fig.  25. 


Target 

ID 

*0  (cm) 

To  (cm) 

Zo  (cm) 

<t>n 

dip  n 

1 

81 -mm 

-23.26 

22.5 

56.16 

0 

-18.3 

2 

105 -mm 

-20.26 

22.5 

69.14 

180 

-18.6 

3 

BLU26 

0.00 

22.5 

43.21 

0 

0.0 

4 

57-mm 

5.22 

22.5 

51.45 

180 

306.6 

5 

60-mm 

0.00 

22.5 

54.50 

0 

0.0 

Table  1 1 :  Correct  answers  for  the  blind-test  data  runs. 


Target 

ID 

Xo  (cm) 

yo  (cm) 

zo  (cm) 

<H°) 

dip  o 

1 

81 -mm 

-21.2 

23.6 

56.56 

0 

-2 

2 

105 -mm 

-6.0 

23.7 

67.57 

170 

-20.8 

3 

BLU26 

0.25 

21.8 

47.55 

*0 

*180 

4 

57-mm 

0.75 

21.6 

58.19 

0 

320.4 

5 

60-mm 

0.67 

19.3 

54.02 

*3 

0 

Table  12:  NSMC  model  inversion  results  for  position  and  orientation  for  blind  test  data.  Numbers 
with  an  asterisk  are  arbitrary  due  to  BOR  considerations. 


A.2  Previous  Multiple  Target  Inversion  Results 

Though  we  would  not  use  the  GEM-3  for  the  research  proposed  here,  our  prior 
work  in  this  area  involves  GEM-3  data. 

We  acquired  data  using  the  GEM-3  instrument  in  2005  over  9  sets  of  targets  in 
the  CRREL  test  plots.  In  each  test  hole,  we  emplaced  either  1  UXO  and  1  piece 
of  clutter,  or  2  UXOs.  We  then  acquired  data  at  two  elevations  over  each  test 
hole  and  used  these  data  in  a  blind  test.  Sample  data  from  one  of  the  test  holes  is 
shown  in  Fig.  26.  As  can  be  seen  from  the  picture,  it  is  not  at  all  clear  by  eye  how 
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Figure  25:  NSMC  model  comparison  to  library  for  blind  test  #3. 


many  targets  are  present  in  the  data  much  less  where  they  are  located.  The  NSMC 
model  was  used  to  simultaneously  invert  for  the  positions,  orientations  and  Q  (or 
the  total  NSMC)  for  each  of  these  test  holes.  Results  are  shown  in  Fig.  27.  The 
NSMC  algorithm  was  able  to  discriminate  between  targets  in  close  proximity  in 
some  cases,  but  in  others,  the  results  were  less  clear. 

While  encouraging  this  study  provides  an  excellent  example  of  the  problems 
faced  when  signals  from  objects  overlap  [11,  14,  68].  These  problems  include 
positional  uncertainties,  scalar  data,  and  sensor  limitations.  The  research  we  pro¬ 
pose  here  aims  to  overcome  these  limitations  with  clutter  mitigation  techniques, 
an  N-targct  locator,  and  more  diverse,  higher  quality  data. 

A.2.a  GEM-3  Data  over  2  UXO 

Data  were  collected  on  a  7x7  grid  at  two  heights  above  each  test  hole.  The  NSMC 
was  then  used  in  an  iterative  fashion  in  an  attempt  to  discriminate  both  UXO  si¬ 
multaneously. 
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ts 

□ 


□stance  £m] 


DEtaice  [ot] 


Figure  26:  GEM-3  data  acquired  over  2  buried  targets. 


A.2.b  Projection  of  GEM-3  Data  over  1  UXO  Obscured  by  Clutter 


Data  projection  involves  using  the  data  from  mono-static  or  multi-static  instru¬ 
ments  to  achieve  the  same  purpose  as  that  for  bistatic  instruments,  that  is,  to  use 
data  at  lower  elevations  (closer  to  clutter)  to  provide  estimates  on  what  data  would 
look  like  if  acquired  at  higher  elevations  (farther  from  clutter).  An  early  attempt 
at  a  data  projection  technique  for  GEM-3  data  is  described  here. 

We  acquired  data  from  a  105mm  UXO  (Fig.  28)  which  was  obscured  by  the 
clutter  layer  affixed  to  a  plywood  substrate  as  shown  in  Fig.  29.  The  clutter  layer 
was  situated  5cm  above  the  upmost  part  of  the  UXO,  and  GEM- 3  measurements 
were  taken  at  2  heights:  15cm  and  20  cm  above  the  upmost  part  of  the  UXO 
respectively.  The  Hz  at  90Hz  from  the  GEM-3  is  shown  at  these  two  height  in 
Fig.  30  (top  subfigures).  Using  our  preliminary  upward  projection  technique,  we 
were  able  to  use  these  data  at  lower  heights  to  estimate  data  as  if  acquired  at 
higher  elevations  (in  this  case  25cm  and  30cm  above  the  UXO).  These  projected 
in  Fig.  30  (bottom  subfigures)  show  the  diminished  influence  of  the  clutter  layer 
as  a  function  of  height.  Note  that  this  clutter  layer  is  an  extreme  case  (to  begin 
with)  and  reacts  almost  as  a  sheet  of  dispersed  clutter,  which  has  a  different  power 
law  than  discrete  targets  [69]. 

The  decay  of  the  projected  data  for  both  the  UXO+clutter  case  and  the  clutter 
only  case  both  approach  the  same  power  law.  This  suggests  that  the  higher  up 
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n 
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Correct 
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Clutter 

Correct 

Correct 

3 

i 

Clutter 

Incorrect 

Correct 

4 

j 

|! 

Clutter 

Incorrect 

Correct 

5 

Correct 

Correct 

J 

|» 

6 

■ 
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a 

1 1 

Correct 

Incorrect 

7 

f 

!* 

|* 

|! 

Incorrect 

\ 

Correct 

8 

An  aluminum 
plate  ?? 

Correct 

Incorrect 

9 

i 

■j 

1 

'1 

Incorrect 

Correct 

Figure  27:  GEM-3  data  inversion  results  from  Partners  poster  2006. 


the  data  is  projected,  the  less  the  influence  the  clutter  has  until  only  the  EMI  re¬ 
sponse  from  the  UXO  is  non-negligible.  This  effect  may  be  exploited  for  use  in 
discrimination  and  needs  further  research. 
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Figure  30:  As  the  sensor  is  raised,  the  contribution  of  clutter  to  the  overall  EMI  response  is  reduced 
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